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Foreword 


The  Fleet  Numerical  Oceanography  Center  is  chartered  to  generate  daily 
forecasts  of  detection  range  for  ASW  sonar  systems  that  operate  in  a  wide 
variety  of  ocean  environments.  This  report  reviews  the  surface-duct  propagation 
loss  prediction  capabilities  of  the  Ship  Helicopter  Acoustic  Range  Prediction 
System. 
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Executive  Summary 


The  propagation  submodels  of  the  Ship  Helicopter  Acoustic  Range 
Prediction  System  (SHARPS)  are  reviewed  within  the  context  of  their 
effectiveness  in  making  accurate  predictions  for  surface-duct  environments. 
The  surface-duct  propagation  loss  calculations  performed  by  SHARPS  are 
based  on  the  submodels  contained  in  Active  RAYMODE.  Fleet  messages  have 
been  critical  of  SHARPS  target  detection  range  predictions  for  in-layer 
and  cross-layer  sonar-target  geometries.  A  review  is  presented  of  those 
computational  methods — common  to  both  Active  and  Passive  RAYMODE— 
most  likely  to  affect  surface-duct  predictions,  and  includes  a  summary  of  the 
main  features  of  a  surface-duct  leakage  subroutine  found  in  one  of  the  previous 
versions  of  Passive  RAYMODE.  A  means  of  coupling  surface-scattered  energy 
into  the  normal  mode  sum  is  reviewed.  The  possibility  of  “anomalous”  surface- 
duct  predictions  is  briefly  addressed,  and  several  ways  to  achieve  improved 
surface-duct  propagation  loss  predictions  are  recommended. 
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Analysis  of  Surface-duct  Propagation  Loss 
Modeling  in  SHARPS 


I.  Introduction 

This  report  presents  a  qualitative  assessment  of  the 
surface-duct  propagation  loss  prediction  capability  of 
the  Ship-Helicopter  Acoustic  Range  Prediction  System 
(SHARPS).  SHARPS  is  used  by  the  Fleet  Numerical 
Oceanography  Center  (FNOC)  to  generate  daily  fore¬ 
casts  of  detection  ranges  for  several  operational  sonar 
systems.  The  Active  RAYMODE  model  generates  the 
basic  acoustic  quantities  used  by  the  latest  version  of 
SHARPS.  The  propagation  loss  subroutines  of  Active 
RAYMODE  were  derived  from  Passive  RAYMODE, 
the  latter  being  the  new  Navy  standard  propagation 
loss  model  for  range-independent  ocean  environ¬ 
ments.  Since  the  installation  of  Active  RAYMODE  in 
SHARPS,  FNOC  has  received  messages  critical  of 
“anomalous”  detection  range  predictions  under  certain 
surface-duct  conditions. 

Passive  RAYMODE  has  been  declared  the  Navy 
standard  propagation  loss  model  for  range-independent 
ocean  environments.  As  a  consequence,  all  propagation 
loss  modules  used  in  generating  fleet  prediction 
products  had  to  be  changed  to  meet  the  new  standard. 
The  decision  was  made  that  even  though  such  changes 
would  be  extremely  costly,  the  costs  were  worth  the 
desired  result:  commonality  of  all  propagation  loss 
modules  used  in  generating  fleet  prediction  products. 
Thus,  both  the  Active  and  the  Passive  versions  of 
RAYMODE  have  been  (or  are  being)  installed  at 
FNOC  to  provide  the  basic  transmission  loss  calcula¬ 
tions  used  in  several  fleet  performance  prediction 
systems. 

The  specific  prediction  system  of  interest  here  is 
SHARPS.  This  system  was  developed  at  FNOC  in  the 
late  1960s.  In  its  original  form,  SHARPS  was  purely 
empirical.  In  the  early  1970s,  SHARPS  was  completely 
revised  in  an  effort  to  incorporate  the  physics  of  under 
water  sound.  An  attempt  to  use  the  original  version 
of  the  Navy  Interim  Surface  Ship  Model  (NISSM) 
failed  due  to  unacceptable  run-times.  A  trimmed 
version  of  NISSM,  referred  to  simply  as  FAST  NISSM 
(Watson  and  McGirr,  1972;  McGirr,  et  al.,  1972),  was 
finally  installed  for  operational  use  in  1973.  This 
code  (then  identified  as  SHARPS  II)  had  many 
shortcomings  and  was  replaced  in  1977  by  a 


“streamlined”  version  of  NISSM  II  (Weinberg,  1973; 
Kirby,  1982).  From  that  point  until  1986  this  version 
of  the  SHARPS  model  (then  identified  as  SHARPS 
III)  underwent  numerous  upgrades. 

One  of  the  major  problems  with  both  SHARPS  II 
and  SHARPS  III  centered  on  the  surface-duct  model 
(based  on  the  AMOS  equations  in  both  versions)  and, 
in  particular,  the  procedure  used  in  defining  the  sonic- 
layer  depth  (for  an  amplification  of  this  problem,  see 
McGirr,  1983).  The  sonic  layer  depth  problem  was 
never  satisfactorily  resolved.  Indeed,  the  analyses  of 
Renner  and  Kirby  (1983)  revealed  that  a  pre-SHARPS 
profile-processing  program  occasionally  yielded 
erroneous  surface-duct  parameters.  This  problem,  in 
combination  with  an  overly  simplified  treatment  of 
surface-duct  propagation,  resulted  in  unacceptably 
large  swings  in  detection  range  whenever  the  near- 
surface  velocity  gradients  were  exceptionally  close  to 
zero. 

Virtually  the  same  problem  exists  with  the  new 
version  of  SHARPS,  except  that  no  special  treatment 
is  given  to  surface-duct  propagation.  The  near-surface 
propagation  loss  prediction  problems  that  existed  in 
previous  versions  of  SHARPS  were  expected  to 
disappear  with  the  installation  of  Active  RAYMODE. 
This  supposition  evidently  emerged  from  expectations 
thatihe  “modal”  character  of  the  new  propagation 
code  would  automatically  take  interduct  coupling  into 
account,  thereby  precluding  any  need  for  a  surface- 
duct  submodel.  These  expectations  were  soon  replaced 
by  expressions  of  concern,  however,  once  the  message 
traffic  revealed  the  reality  of  the  situation:  the  “normal 
mode”  method  used  in  Active  RAYMODE  essentially 
ignores  the  coupling  of  energy  from  one  duct  to 
another.  Thus,  the  surface-duct  problem  persists,  and 
FNOC  personnel  once  again  find  themselves  trying  to 
respond  to  criticisms  from  the  Fleet  regarding  the 
“same  old  problem.” 

Essentially,  recipients  of  SHARPS  forecasts  question 
in-layer  and  cross-layer  detection  range  predictions.  For 
presumably  well-defined  surface  ducts,  predicted 
detection  ranges  for  both  sonar  and  target  in  the  duct 
are  shorter  than  detection  ranges  predicted  for  sonar 
in  the  duct  and  target  below  the  duct.  However, 
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on-board  predictions  for  the  same  inputs  yield  the 
opposite  results.  The  first  question  is,  were  identical 
inputs  used  by  both  the  ashore  and  the  afloat  prediction 
systems?  The  second  question  is,  are  both  the  ashore 
and  the  afloat  implementations  of  the  basic  acoustic 
model  (namely,  Active  RAYMODE)  identical? 

The  first  question  can  be  answered  only  by  Fleet 
and  FNOC  personnel  responsible  for  making  the 
predictions.  The  second  question,  however,  brings 
up  implementation  issues.  Both  Active  and  Passive 
RAYMODE  are  under  configuration  management  and, 
as  part  of  each  software  package,  test  cases  are  included 
for  verification.  Presuming  that  the  Active 
RAYMODE  software  was  verified  after  being  hard¬ 
wired  into  the  SHARPS  “shell”  and  that  a  similar 
verification  procedure  is  faithfully  followed  after 
installation  in  each  system  afloat,  the  logical  course 
of  action  entails  comparing  surface-duct  predictions 
generated  by  ashore  and  afloat  version  of  Active 
RAYMODE  against  high-confidence  control  models. 
If  acceptable  agreement  is  not  obtained  by  either 
version,  then  the  only  recourse  is  to  analyze  the  physics 
contained  in  Active  RAYMODE  that  deals  with 
surface-duct  environments. 

The  verification  procedure  alluded  to  above  may 
not  uncover  “bugs”  peculiar  to  implementations  on 
specific  machines  or  operating  systems.  Steps  have 
been  taken  by  both  the  developing  agency,  the  Naval 
Underwater  Systems  Center  (NUSQ,  and  the  organiza¬ 
tion  responsible  for  configuration  management.  Naval 
Oceanographic  Office  (NAVOCEANO),  to  ensure  that 
there  is  no  machine/operating-system  dependency.  The 
RAYMODE  codes  delivered  to  the  Naval  Ocean 
Research  and  Development  Activity  (NORDA)  were 
programmed  in  strict  FORTRAN  77,  which  should  be 
free  of  any  such  dependencies.  Nonetheless,  there  is 
always  the  possibility  that  a  model  which  generates 
acceptable  results  on  the  HP-9020/UNIX  will  not, 
for  some  cases,  generate  identical  results  on  a  VAX 
running  VMS.  Even  though  there  may  be  some  ques¬ 
tions  regarding  verification,  which  are  presumably 
being  addressed  by  the  various  software  firms  involved, 
the  purpose  of  the  effort  reported  herein  is  to  uncover 
weaknesses  in  the  Active-RAYMODE  model  with 
regard  to  surface-duct  predictions. 

The  report  is  broken  down  into  four  major  sections. 
Section  II  reviews  those  aspects  of  the  RAYMODE 
physics  modeling  that  might  have  some  bearing  on 
surface-duct  predictions.  Section  III  addresses  surface- 
duct  propagation  modeling  problems  in  general. 
Section  IV  examines  the  impact  of  neglecting  surface 
scattered  energy.  Section  V  illustrates  jump  discon¬ 
tinuities  that  are  common  to  methods  based  on  hybrid 
solutions,  and  Section  VI  closes  the  report  with 
summary  remarks  and  recommendations. 


II.  A  Review  of  RAYMODE  Model 
Physics 

Passive  RAYMODE  was  initially  conceived  in  the 
late  1960s  by  Leibiger  (1968)  and  has  evolved  into  a 
widely  used  propagation  code.  In  essence,  Leibiger 
replaces  the  single  contour  integral  of  the  Fourier- 
Bessel  solution  to  the  reduced  wave  equation  into 
several  integrals,  each  of  which  can  be  associated  with 
a  family  of  rays.  These  integrals,  which  result  from 
expanding  the  reciprocal  of  the  Wronskian  into  an 
infinite  series,  are  then  evaluated  using  an  approximate 
normal-mode  method  if  the  number  of  modes  is  not 
too  large  or  otherwise  using  integration  techniques 
more-or-less  unique  to  RAYMODE.  The  advantages 
offered  by  RAYMODE  over  standard  normal-mode 
treatments  are  two-fold.  First,  the  computational 
load  does  not  significantly  increase  with  increasing 
frequency,  as  is  the  case  with  normal-mode  models. 
Second,  the  solution  is  partitioned  into  components 
that  give  rise  to  plausible  geometrical  interpretations 
similar  to  those  available  from  methods  based  strictly 
on  ray  acoustics. 

The  propagation-loss  algorithms  used  in  Active 
RAYMODE  essentially  form  a  subset  of  those  used 
in  Passive  RAYMODE.  As  a  consequence.  Active 
RAYMODE  users  are  concerned  about  the  loss 
in  accuracy  that  could  result  from  calculations 
restricted  to  the  high-frequency  algorithms  of  Passive 
RAYMODE.  Concerns  of  the  same  sort  also  pertain 
to  Passive  RAYMODE,  since  numerous  assumptions 
and  approximations  have  been  made  in  an  effort  to 
obtain  fast  execution  times.  The  focus  here  is  on  those 
aspects  of  RAYMODE  that  potentially  impact  surface- 
duct  predictions. 

A.  Methods  Used  in  Passive  RAYMODE 

Certain  documents  (Leibiger,  1968,  1971; 
Deavenport,  1978;  Davis  and  Council,  1985a,  1985b) 
present  various  aspects  of  the  model  physics.  The  text 
by  DiNapoli  and  Deavenport  (1979)  is  an  excellent 
reference  on  contemporary  propagation  modeling  in 
general,  and  also  presents  a  brief  discussion  of 
RAYMODE.  Those  documents  written  by  Davis  and 
Council  are  the  most  complete,  although  the  sections 
that  address  the  special  surface-duct  treatment  do  not 
apply  to  the  most  recent  version  of  Passive 
RAYMODE,  since  the  surface-duct  module  is  not 
included  in  the  latest  version.  Moreover,  this  surface- 
duct  treatment  evidently  has  never  been  considered  for 
inclusion  in  Active  RAYMODE. 

The  steps  leading  to  the  RAYMODE  versions  of 
normal-mode  and  multipath  expansion  solutions  are 
developed  in  Appendix  A.  Much  of  the  material 
presented  in  Appendix  A  parallels  the  Passive 
RAYMODE  documentation  of  Davis  and  Council 
(1985b).  Readers  who  are  interested  in  understanding 
the  model  physics  should  peruse  all  of  the  references. 


For  surface-duct  paths,  three  RAYMODE  methods 
can  be  applied  to  solve  the  multipath  expansion 
integral,  although  the  surface-duct  option  offered  in  a 
previous  version  of  RAYMODE  essentially  represents 
a  fourth  RAYMODE  method.  The  other  three  methods 
are  (1)  numerical  integration  based  on  the  fast,  Fourier 
transform  (FFT);  (2)  the  original  RAYMODE  inte¬ 
gration  (ORI)  method;  and  (3)  a  special  high-frequency 
integration  (HFI)  method.  If  the  user  does  not  select  a 
particular  method,  then  the  program  selects  one. 
For  frequencies  less  than  3000  Hz,  either  of  the  first 
two  methods  is  selected,  depending  on  estimated 
execution  time  for  each  method.  If  the  best  option  is 
the  FFT  method,  but  aliasing  is  anticipated,  then  the 
ORI  method  is  used  instead.  For  frequencies  greater 
than  3000  Hz,  the  special  HFI  method  is  selected.  The 
methods  of  primary  interest  here  are  the  original 
RAYMODE  integration  method  and  the  special  HFI 
method.  The  special  surface-duct  treatment  that 
Leibiger  applied  in  a  previous  version  of  RAYMODE 
is  discussed  in  Section  III.  The  ORI  and  HFI  methods 
are  discussed  in  the  next  two  subsections. 

B.  Original  RAYMODE  Integration 

The  general  form  of  the  RAYMODE  multipath 
expansion*  integral  (see  Eq.  (A38)  in  Appendix  A)  is 

/(*„,*„)  =JW)  exp {-i[kr  +  a(*)l}  dk,  (1) 

where  the  limits  of  integration  extend  from  ka  to  kb. 
An  approach  often  used  in  the  evaluation  of  integrals 
of  this  form  is  the  method  of  stationary  phase.  This 
method  is  summarized  in  Appendix  B. 

Leibiger  interrupts  the  stationary-phase  analytical 
process  and  does  not  actually  proceed  to  the  final 
formula  (see  Eq.  (B51)  of  Appendix  B).  The  reason 
for  this  circumvention  is  that  the  conditions  required 
to  achieve  the  final  step  may  not,  for  all  frequencies 
of  interest,  be  met.  An  acceptable  application  of  the 
formula  requires  that  the  function  in  the  exponent, 
h(k),  be  multiplied  by  a  large  constant.  The  constant 
in  the  development  as  presented  in  Appendix  B  is  unity. 
The  argument  of  the  phase  function  can  be  expressed 
as,  for  example,  k0[kr/k0  +  a(k)/k0\,  and  the 
variable  of  integration  transformed  to  k/k0,  where  k0 
(=  2 nf/C0)  is  some  reference  wavenumber.  Thus, 
the  stated  condition  is  met  for  high  frequencies. 
For  low  frequencies,  the  limits  of  integration  must 
be  confined  to  a  relatively  small  interval,  where 
(hopefully)  the  resulting  integral  yields  an  accurate 
estimate.  Early  work  by  Leibiger  (1968)  indicates  that 
the  stationary-phase  formula  was  initially  exploited  and 


•For  an  elementary  interpretation  of  this  technique  consult 
Sect.  35  in  the  text  by  Brekhovskikh  (1980). 


then  discarded  in  favor  of  the  method  outlined  in 
Appendix  B. 

Leibiger  (1971)  demonstrates  the  relationship 
between  the  stationary-phase  “formula”  and  the 
standard  ray-acoustical  expression.  In  an  earlier  report 
(Leibiger,  1968),  he  demonstrates  that  the  existence  of 
a  point  of  stationary  phase  implies  the  existence  of  a 
ray  trajectory  connecting  the  source  and  the  receiver. 

The  ORI  method,  as  presented  in  Appendix  B, 
suggests  the  not  entirely  implausible  possibility  that  the 
amplitude  and  phase  factors  do  not  necessarily  accede 
to  conditions  generally  assumed  extant  in  justifying  the 
application  of  stationary-phase  methods.  These  matters 
are  discussed  in  Appendix  C.  Even  when  the  evaluation 
occurs  at  an  end  point,  a  stationary  phase  analysis  can 
be  used  successfully  to  determine  asymptotic  behavior. 
A  precise  assessment  of  the  impact  of  asymptotic 
behavior  is  beyond  the  scope  of  this  effort,  but  as  a 
practical  matter,  in  any  problem  where  a  function 
is  to  be  approximated,  some  a  priori  upper  bound  is 
placed  on  the  tolerable  error.  A  perusal  of  the  computer 
code  indicates  that  Leibiger  has  taken  precautionary 
steps  to  circumvent  serious  problems. 

Some  hint  about  the  magnitude  of  error  attributable 
to  the  ORI  method  can  be  gleaned  from  Bartberger’s 
(1981)  evaluation  of  a  previous  version  of  Passive 
RAYMODE.  Bartberger  notes  that  differences  between 
the  normal-mode  method  and  the  ORI  method  (which 
may  have  been  altered  since)  can  be  significant  for 
coherently  summed  outputs  but  appear  to  be  minor 
for  incoherently  summed  outputs.  Discontinuous 
jumps  that  are  created  when  the  method  of  solution 
switches  from  one  form  to  another  are  briefly  examined 
in  Section  IV. 

Leibiger  (1971)  notes  that  his  particular  stationary- 
phase  treatment  avoids  certain  problems  that  can 
develop  as  the  acoustic  frequency  increases.  At  a  high 
enough  frequency,  the  ORI  method  is  supplemented 
by  a  special  high  frequency  integration  HFI  approach. 

C.  Special  High-frequency  Integration 

Keep  in  mind  that  a  derivation  of  the  high-frequency 
procedure  has  not  been  documented  by  the  model 
developer,  so  the  only  information  available  is  the  brief 
description  given  by  Davis  and  Council  (1985b),  and 
the  computer  code.  The  main  features  of  this  technique 
are  presented  in  Appendix  B. 

A  perusal  of  the  computer  code  for  details  on  how 
the  HFI  procedure  is  implemented  reveals  that  the  Airy 
integral  is  approximated  by  a  polynomial  of  degree  six 
when  the  argument  has  magnitude  less  than  five,  and 
by  the  inverse  square  of  the  argument  otherwise.  This 
particular  form  of  asymptotic  representation  seems 
rather  crude,  although  quite  possibly  justified.  The 
polynomial  approximation  appears  to  be  a  truncated 
power  series  expansion,  which  probably  can  be 


improved  upon  by  using  an  economized  rational 
polynomial  approximation  instead. 

D.  Surface  Reflection  Loss 

RAYMODE  surface  reflection  loss  calculations  are 
based  on  an  ad  hoc  adaptation  of  results  derived  in 
part  from  theoretical  considerations  and  in  part  from 
experimental  data  (Marsh  and  Schulkin,  1962).  A 
summary  of  this  submodel  is  presented  in  Appendix  D. 

III.  Surface-duct  Propagation  Modeling 

This  section  addresses  several  more  or  less 
independent  aspects  of  modeling  propagation  loss  for 
surface-duct  environments.  One  of  the  working 
hypotheses  upon  which  this  report  is  based  is  that  the 
surface-duct  model  used  in  the  1983  version  of 
RAYMODE  can,  with  suitable  modifications,  yield 
improved  surface-duct  predictions  if  incorporated  into 
the  latest  versions  of  the  RAYMODE  codes.  Modifica¬ 
tions  are  necessary  for  two  reasons:  there  appear  to 
be  errors  in  the  computer  code  and  some  means 
should  be  incorporated  to  account  for  the  coupling  of 
surface-scattered  energy  into  the  below-layer  region. 
This  latter  topic  is  addressed  in  Section  IV. 

The  RAYMODE  special  surface-duct  algorithms 
alluded  to  are  discussed  in  Section  A.  Section  B  presents 
some  characteristics  of  the  two-segment  n2-linear 
velocity  profile  model,  which  forms  the  basis  of  the 
surface-duct  module  used  in  the  latest  versions  of 
FACT  and  ASTRAL.  This  approach  to  the  problem 
has  both  advantages  and  disadvantages.  A  procedure 
that  was  successfully  used  in  the  1960s  is  reviewed,  and 
several  detracting  features  of  the  /i2-linear  model  are 
discussed. 

A.  RAYMODE  Surface-duct  Model 

The  normal-mode  sum  used  in  Active  RAYMODE 
is  possibly  adequate  to  handle  range-invariant,  shallow- 
water  ducts.  For  deep-ocean  profiles  that  contain 
multiple  ducts,  there  is  no  provision  for  cross-channel 
coupling.  Also,  the  mode  stun  is  limited  to  a  maximum 
of  10  propagating  modes,  so  for  frequencies  above 
about  3  kHz  most  ducts  are  likely  to  be  handled  by 
the  HFI  method.  In  many  normal-mode  codes,  the 
modes  (eigenvalues)  are  determined  by  numerical 
integration  over  the  entire  water  column.  Some  codes 
use  iterative  procedures  to  solve  the  characteristic 
equation.  In  either  case,  interchannel  coupling  effects 
are  automatically  included.  The  simplified  procedure 
used  in  Active  RAYMODE,  as  it  stands,  cannot 
address  this  type  of  coupling. 

The  1983  version  of  Passive  RAYMODE  includes 
the  option  of  a  special  treatment  of  the  normal-mode 
method  when  applied  to  surface  ducts.  In  the  standard 
RAYMODE  approach,  the  magnitude  of  the 
reflection  coefficient  at  each  k-segment  interface  is  one. 


except  at  the  surface  and  the  bottom.  In  the  special 
treatment  accorded  surface  ducts,  generalized  Wentzel- 
Kramers-Brillouin  (WKB)  reflection  coefficients  are 
introduced.  These  reflection  coefficients  are  then  used 
in  constructing  complex  eigenvalues,  where  each 
imaginary  component  corresponds  to  an  exponential 
loss  of  energy  with  range.  Thus,  for  a  source  located 
within  the  duct,  WKB  modes  lose  some  of  their  energy 
at  each  turning  point,  where  the  lost  energy  is 
transmitted  into  the  region  below  the  layer  depth. 
Similarly,  modes  excited  by  a  source  located  in  the 
negative-gradient  region  transfer  some  of  their  energy 
to  waves  that  propagate  within  the  duct.  This 
procedure,  if  added  to  the  present  version  of  Active 
RAYMODE,  should  yield  improved  estimates  of 
surface-duct  propagation.  To  have  any  impact  at 
frequencies  above,  say,  3  kHz,  the  maximum  number 
of  modes  would  have  to  be  extended  to  about  50.  Thus, 
questions  naturally  arise  pertaining  to  run-time  and 
accuracy  required  in  the  determination  of  the  complex 
eigenvalues. 

The  only  documentation  on  the  basis  for  this  special 
surface-duct  procedure  is  the  brief  account  presented 
by  Davis  and  Council  (1985b).  (The  generalized  WKB 
reflection  coefficient  procedure  is  described  in  the 
article  by  Murphy  and  Davis,  1974)  The  modifications 
to  the  mode  sum  depend  on  source-receiver  geometry. 
When  both  the  source  depth  (z,)  and  the  receiver 
depth  (zr)  are  less  than  the  layer  depth  (zL),  the  mode 
sums  are  applicable  using  appropriate  reflection 
coefficients  at  the  turning  points.  The  upper  turning 
point  in  this  case  is  the  surface;  therefore,  Ru 
corresponds  to  the  surface  reflection  coefficient 
(discussed  in  Appendix  A).  At  the  lower  turning  point 
depth,  say,  z,,  R,  is  replaced  by  RA^z,),  a  generalized 
WKB  reflection  coefficint.  The  subscript,  AB. 
indicates  the  sense  of  direction  associated  with  the 
coefficient,  where  A  stands  for  above-layer  and  B 
stands  for  below-layer.  Details  pertaining  to  the 
calculation  of  this  coefficient,  along  with  the  corres¬ 
ponding  transmission  coefficient,  are  presented  in 
subsequent  paragraphs. 

When  both  the  source  and  the  receiver  are  situated 
below  the  layer  depth,  the  reflection  coefficient  at  the 
upper  turning  point  is  replaced  by  /^(zj.  which  is 
similar  to  RAB(z,)  for  the  in-layer  case,  but  where 
the  above-layer  and  below-layer  parameters  are 
interchanged.  The  lower  turning  point  is  assumed  to 
be  at  depth;  therefore,  R,  is  assumed  to  have  a 
magnitude  of  1  and  a  phase  of  -n/2. 

For  the  cross-layer  case,  say,  zs  <  zL  and  zr  >  zL, 
Leibiger  takes  a  different  approach.  Instead  of 
reflection  coet'ficients,  the  amplitude  factors  are 
multiplied  by  a  transmission  coefficient.  Thus,  energy 
that  would  otherwise  be  considered  trapped  using  WKB 
modes,  if  indeed  they  even  existed,  now  would  have 
the  potential  of  leaking  into  the  below-layer  region. 


The  “mode-leakage”  phenomenon  is  perhaps  best 
explained  by  drawing  upon  ray  analogies.  Figure  1 
illustrates  a  ray  trajectory  that  turns  around  at  the  layer 
depth.  To  indicate  that  some  of  the  energy  carried  by 
the  ray  remains  within  the  duct  and  that  some  leaks 
into  the  negative  gradient  region  below,  the  ray  is 
split  into  two  components  at  the  turning  point. 

To  extend  the  ray  analogy  a  bit  further,  there  are 
partial  reflections  at  the  ray  turning  points.  If  Ut 
and  Vt  denote,  for  layer  two  linearly  independent 
solutions  of  the  depth-separated  wave  equation,  then 
at  the  layer  depth  (see  Bucker,  1980), 

£/,  +  R  y,  =  T  U2,  (2) 

where  (/,  and  Vx  denote  down-going  and  up-going 
waves  within  the  duct  (layer  1)  and  U2  denotes  a 
down-going  wave  in  the  region  below  the  duct 
(layer  2).  The  reflection  and  transmission  coefficients 
at  the  bottom  of  the  duct  may  be  obtained  from  Ryan 
(1980): 

|/?|2  =  (1  +  s3)/(l  +  s)3  (3) 

and 

i-i*i*,  (4) 


Figure  I.  Ray  trajectory  splitting  at  layer  depth. 


where  s3  =  P/f3 2.  The  parameters  /?,  and  ft2  are 
the  gradients  of  the  squared  index  of  refraction 
for  each  layer  (see  Fig.  2).  Since  1  +  si  =  (1  +  s) 
(1  -  5  +  s2),  the  expression  for  \R\  can  be  written 

*  =  (1  -  s  +  s2\ ~/( 1  +  s).  (5) 

Thus,  the  transmission  coefficient  is  given  by 

|T|  =  (3s);V(l  +  s).  (6) 

These  expressions  for  \R\  and  \T\  are  in  the  forms 
presented  in  the  RAYMODE  documentation  of  Davis 
and  Council  (1985b).  Coefficients  of  the  same  basic 
form  can  be  found  in  the  1983  RAYMODE  computer 
code,  although  the  definition  of  s  found  there  is 
proportional  to  the  reciprocal  of  the  one  defined  here.* 

Figure  3  indicates  the  ways  (excluding  surface 
scattered  paths)  in  which  energy  can  propagate  from 
an  in-layer  source  to  a  below-layer  receiver.  Figure  3 

(a)  depicts,  in  terms  of  ray  equivalents,  an  in-layer 
trapped  mode  and  a  below-layer  untrapped  mode. 
These  modes  interact  by  means  of  diffraction  leakage 
through  the  barrier  of  thickness' d  =  z,-zu.  Figure  3 

(b)  illustrates  the  split-beam  ray,  corresponding  to 
transitional  modes,  or  those  modes  that  are  neither 
strongly  trapped  nor  decidedly  leaky.  Figure  3  (c)  shows 
a  steep  ray  that,  in  standard  ray-acoustical  treatments, 
transits  directly  into  the  below-layer  region. 

The  inclusion  of  a  reflected  component  is  indicative 
of  the  situation  associated  with  a  leaky  mode.  Davis 
and  Council  (1985v)  point  out  that  a  generalized 


•This  circumstance  is  a  good  example  of  why  the  physics  of 
a  model,  especially  one  intended  for  operational  applications, 
should  be  fully  documented  and  subjected  to  peer  review. 


Figure  2.  An  n2-linear  profile  for  two-layer  surface  duct  (a)  sound-velocity  profile  and  (b)  mdex-of-refraction 
profile. 
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WKB  phase  integral  approach,  if  implemented 
properly,  can  be  used  to  simulate  the  interduct  energy- 
coupling  process.  The  reader  interested  in  more  details 
concerning  generalized  WKB  procedures  should  consult 
Murphy  and  Davis  (1974).  The  approach  used  in 
Passive  RAYMODE  is  somewhat  involved,  but 
essentially  a  curve  is  developed  that  gives  \R\,  and 
therefore  |  T\ ,  as  a  function  of  in-layer  and  below-layer 
parameters  defined  by 


discourage  all  but  the  most  persistent  from  pursuing 
this  course  of  action. 

SDUCT  is  highly  convoluted  and  difficult  to 
interpret.  Unraveling  parts  of  this  code  reveals  some 
curious  expressions.  For  example,  one  line  of  the 
SDUCT  subroutine  reads 

RL  =  1.  -  HLD  +  HLD*RLP. 


(in-layer)  x,  =  h](k2L  -  k2J, 


Tracing  one  particular  path  through  the  code  leads  to 
the  following  interpretation  of  this  statement: 


and 

(below-layer)  x2  =  h22(k2  -  Arp*  (8) 

where  km  =  c o/Cm,  kL  =  o o/CL,  A,  =  (Ar^/?,)-*,  and 

h2  =  (k2J2y'A. 

For  untrapped  modes  (Fig.  3  (c);  Ryan,  1980)  the 
reflection  coefficient  can  be  expressed  as 

1*1  =  1(1  +  s3)/(4*,m  i  *  1,2.  (9) 

A  curve-fit  can  then  be  made  using  points  generated 
by  this  equation  evaluated  at  values  of  km  far 
removed  from  ko  at  the  point  determined  for  trapped 
modes,  and  at  a  point  (0  <  z  <  *u)  within  the  duct 
using  results  based  on  a  WKB  approximation. 

Specific  details  of  this  surface-duct  submodel  are  not 
available  in  a  “technical  description”  document  as  is 
customary,  and  hence  the  curious  model  user  must 
resort  to  snooping  through  the  computer  program 
source  code.  Ordinarily,  th  procedure  leads  to 
satisfactory  explanations  corurning  how  the  model 
developer  surmounts  various  implementation 
problems.  In  the  1983  version  of  Passive  RAYMODE, 
the  surface-duct  calculations  are  done  in  subroutine 
SDUCT.  A  brief  review  of  this  code,  however,  will 


RL  =  1  -(/?, -  l)exp j -4/3  2-(l+-,) 
L  Si  P3 


k 

r-1 

*i 


3/2 


(10) 

where  R,  -  [1  -  p  +  p2]V(  1  +  p),  p3  =  -g2/gx,  g, 
and  g2  are  the  in-layer  and  below-layer  velocity 
gradients,  A^,is  the  wave  number  for  mode  m,  and  k, 
is  the  wave  number  associated  with  the  layer  depth. 
The  expression  for  R,  is  identical  in  form  to  the 
expression  given  previously  for  trapped  modes  (except 
that  p  *  -1  /s). 

Using  the  C-linear  velocity  profile  model,  the  last 
factor  in  the  exponent  can  be  interpreted  as 


(f,/CJ*Gi  -  O'" 


(ID 


where  z,  is  the  layer  depth  and  (zm,Cm)  is  the  depth- 
velocity  pair  at  the  turning  point  for  mode  m.  When 
(g,/Cm)3/2  is  combined  with  the  leading  factors  in  the 
exponent,  this  expression  can  be  put  into  the  form 

*£  =  !-(/?,-  1) exp {-V,km p«(  1  -  s3)(z,  -  zj3 , 

(12) 
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where  (3m  =  2gt/Cm.  The  exponent  in  this  last  expres¬ 
sion  bears  a  curious  resemblance  to  exponents  typical 
of  WKB  phase-integral  functions,  except  that  k0f. ?,* 
is  usually  found  in  place  of  kmjim/l.  The  dependence 
of  this  factor  on  mode  number  is  unexpected  and 
deserves  explanation. 

Although  some  of  the  WKB-type  expressions  found 
in  subroutine  SDUCT  are  similar  in  basic  form  to 
those  found  in  published  work  (e.g.,  Barnard 
and  Deavenport,  1978;  Hall,  1976;  Kibblewhite  and 
Denham,  1965),  there  are  enough  differences  to 
warrant  separate  documentation.  The  addition  of  an 
upgraded  version  of  SDUCT  needs  to  be  considered, 
because  the  inclusion  of  generalized  WKB  reflection 
and  transmission  coefficients  is  necessary  to  obtain 
realistic  surface-duct  predictions  using  the  RAYMODE 
normal-mode  sum  (but,  of  course,  with  more  modes). 

B.  Characteristics  of  the  Bilinear  Profile  Model 

The  special  surface-duct  treatment  reviewed  in 
Section  A  makes  use  of  the  nMinear  velocity  profile 
model  to  obtain  simple  expressions  for  the  generalized 
WKB  reflection  and  transmission  coefficients.  This 
profile  model  has  enjoyed  wide  use  in  propagation 
codes,  including  those  based  on  ray  acoustics  and 
normal-mode  theory.  Among  normal-mode  codes  that 
rely  on  analytical  solutions  to  the  z-separated  wave 
equation  for  each  profile  layer,  the  /^-linear  model 
yields  closed-form  solutions  that  can  be  expressed  in 
terms  of  Airy  functions  or  modified  Hankel  functions 
of  order  one-third. 

This  profile  model  also  results  in  straightforward 
evaluations  of  range  and  travel-time  integrals. 
Although  the  C-linear  (constant  gradient)  profile  model 
accedes  to  simple  ray-acoustical  expressions,  it  does 
not  lend  itself  quite  so  readily  to  wave-acoustical  com¬ 
putations.  Essentially,  the  nMinear  profile  model 
represents  a  mathematically  expedient  approach  to 
solving  wave-acoustical  problems,  and  is  probably 
satisfactory  for  many  modeling  requirements.  Some 
of  the  advantages  and  disadvantages  of  using  models 
based  on  closed-form  solutions  vice  numerical  integra¬ 
tion  are  reviewed  in  subsection  B.l. 

Two  features  of  the  /^-linear  profile  model  detract 
from  its  appeal,  at  least  in  regard  to  its  appropriateness 
for  surface  ducts.  One  feature  is  the  direction  of 
curvature,  which  is  sometimes  just  the  opposite  of  that 
suggested  by  measured  profile  data;  the  other  is  the 
extreme  discontinuity  artificially  introduced  at  the  layer 
depth.  This  last  feature  is  seldom  evident  in  actual 
profiles.  The  impact  of  these  features  is  discussed  in 
subsections  B.2.-B.4. 

1.  Computational  Pros  and  Cons 

The  bilinear  normal-mode  model  was  initially 
developed  by  Furry  (1945);  also  see  Freehafer,  1951) 
to  solve  tropospheric  electromagnetic  wave  propagation 


problems.  It  was  adapted  for  underwater  acoustics  by 
Marsh  (1950)  and  later  refined  by  Pedersen  and 
Gordon  (1965).  The  Pedersen-Gordon  model  was 
subsequently  exploited  by  Watson  and  McGirr  ( 1966) 
as  an  integral  component  of  a  production-oriented 
performance  prediction  capability.  The  primary  cus¬ 
tomers  were  local  (U.S.  Navy  Electronic  Laboratory— 
NEL)  developers  of  active  sonar  systems.  At  that  time 
computer  technology  had  not  advanced  to  a  stage 
where  numerical  integration  of  the  z-separated  wave 
equation  was  considered  a  viable  option.  The  two-layer 
normal-mode  model,  however,  offered  some  promise 
as  a  computationally  practicable  approach  to  solving 
the  surface-duct  propagation  problem. 

The  characteristic  equation  associated  with  this 
model  takes  the  form  of  a  3x3  matrix  having  elements 
expressible  in  terms  of  modified  Hankel  functions 
of  order  one-third*.  These  functions  were  the  object  of 
intensive  investigation  by  personnel  at  the  Harvard 
University  Computation  Laboratory  (1945)  during 
World  War  II.  Their  efforts  produced  extensive  tables 
valuable  for  verifying  independently  developed 
computer  algorithms,  which  provide  useful  informa¬ 
tion  regarding  power-series  and  asymptotic-series 
representations  of  the  Hankel  functions.  Thus,  with 
much  of  the  computational  ground  work  already 
accomplished,  incentive  to  develop  a  semiautomated 
surface-duct  propagation  modeling  capability  was 
strong. 

The  Pedersen-Gordon  model  (also  referred  to  as  the 
zero-limit  profile  model;  Hall,  1982)  requires  that  two 
families  of  modes  be  considered.  The  second  family 
of  modes,  however,  becomes  an  important  considera¬ 
tion  only  at  short  ranges,  at  low  frequencies,  and  under 
weak  gradient  conditions.  In  this  sense,  the  second 
family  plays  a  role  similar  to  that  of  the  branch-line 
integral  associated  with  other  branch-cut  approaches. 
For  most  applications,  only  the  first  family  of  modes 
is  considered.  The  procedure  adopted  at  NEL  during 
the  late  1960s  entailed  precalculating  first-family  eigen¬ 
values  for  selected  sets  of  the  controlling  parameters. 
The  controlling  parameters  are  p  =  (gA/gB)'’,  where 
gA  and  gB  are  the  above-layer  and  below-layer 
velocity  gradients  evaluated  at  the  top  of  each  layer, 
and  Af,  the  so-called  “ducting”  parameter,  defined  as 

M  =  2n¥'gAv'fv'ZL/Cs .  (13) 

Thus,  for  a  given  value  of  p,  the  real  and  imaginary 
parts  of  the  eigenvalues  were  calculated  and  stored  for 
values  of  M  ranging  from  1  to  100  in  steps  of  0.25. 
Needless  to  say,  even  for  only  60  modes,  this  procedure 
resulted  in  the  storage  of  thousands  of  punched  cards. 


•In  contemporary  literature.  Airy  functions  are  usually  used 
instead. 
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Sets  of  eigenvalues  were  calculated  and  stored  for  only 
a  few  values  of  p,  and  predictions  had  to  be  limited 
to  above-layer  and  below-tayer  gradients  that  resulted 
in  values  of  p  very  close  to  the  few  stored  values.  The 
procedure  for  making  a  prediction  entailed  selecting 
st  -d  sets  of  eigenvalues  for  values  of  M  on  either 
sic.  jf  that  value  determined  from  input  data.  These 
stored  sets  were  then  interpolated  by  the  transmission 
loss  program  to  arrive  at  the  correct  set  of  eigenvalues. 
Since  the  eigenvalues  are  independent  of  source-receiver 
geometry,  predictions  could  be  made  for  several  source- 
receiver  depth  combinations  without  having  to 
recalculate  new  sets  of  eigenvalues.  This  procedure  was 
reasonably  efficient  when  predictions  were  made  for 
a  particular  sonar  system,  i.e.,  a  single  frequency. 
When  predictions  had  to  be  made  at  several  frequen¬ 
cies,  the  procedure  became  labor  intensive. 

This  approach  may  seem  crude  in  comparison  to 
contemporary  prediction  procedures.  Indeed,  such 
models  as  ASTRAL  and  RAYMODE  perform  all 
required  computations  “on  the  fly,”  although  many 
liberties  have  been  taken  by  the  model  developers  in 
making  approximations  to  yield  relatively  fast  execu¬ 
tion  times.  Kuperman  et  al.  (1986)  and  Porter  et  al. 
(1987)  are  developing  the  Wide-area  Rapid  Acoustic 
Prediction  (WRAP)  system  in  which  certain  modal 
quantities  are  precalculated  at  selected  spatial  nodes 
in  an  effort  to  produce  a  high-speed  three-dimensional 
sound  field  prediction  capability.  Thus  the  idea  of 
precalculating  intermediate  modal  quantities  is  still 
considered  a  viable  approach  in  circumventing  calcula¬ 
tions  that  are  numerically  intensive. 

There  is  strong  temptation  to  suggest  that  a  normal- 
mode  model,  in  conjunction  with  precalculated  sets  of 
eigenvalues,  be  considered  for  implementation  at 
FNOC  to  handle  the  surface-duct  propagation  predic¬ 
tion  function.  The  large  storage  capacity  needed  to 
accommodate  thousands  of  precalculated  eigenvalues 
would  not  be  a  serious  imposition  for  shore-based 
mainframes.  The  propagation  modeling  commonality 
requirement  recently  issued  would  necessitate  that  the 
same  prediction  function  be  implemented  in  on-board 
systems  as  well.  Computer  hardware  presently  being 
used  by  on-board  systems  precludes  this  possibility. 
The  addition  of  an  optical  disk  would  lend  credence 
to  a  precalculated  approach,  even  for  on-board 
systems. 

This  approach  has  considerable  appeal  from  a 
strictly  computational  point  of  view.  Some  disadvan¬ 
tages,  however,  need  to  be  considered.  One  of  these 
disadvantages  pertains  to  the  possible  inflexibility  of 
a  two-layer  model.  Experiments;  aata  (Anderson  and 
Pedersen,  1976)  indicate  that  other  distinct  near-surface 
velocity  profile  configurations  are  just  as  prominent 
as  the  surface  duct.  The  modeling  requirements 
imposed  by  these  other  profile  configurations 
demand  a  more  rigorous  multiple-layer  treatment. 
Rigorous  solutions  are  not  likely  to  b^implemented 


without  problems.  The  reader  interested  in  further 
details  on  problems  that  can  be  encountered  is  referred 
to  Pedersen  and  McGirr  (1982,  pp.  97-108). 

Even  when  a  two-layer  modeling  approach  is 
applicable,  there  are  some  concerns  associated  with  the 
n2-linear  profile  model.  These  concerns  are  discussed 
in  the  next  few  sections. 

2.  Sharp  Discontinuity  at  the  Layer  Depth 

The  artificiality  of  the  sharp  discontinuity  that  the 
nMinear  profile  model  imposes  at  the  layer  depth  has 
been  a  source  of  some  concern  (Furry,  1946)  since  the 
model  was  first  exploited.  In  many  instances  actual 
velocity  profile  data  exhibit  what  appears  to  be  a  sharp 
discontinuity.  In  just  as  many  if  not  more  cases,  the 
reversal  in  sign  of  the  gradient  is  gradual.  The  concept 
of  layer  depth  is  indeed  elusive.  Furthermore,  even 
when  a  point  can  be  identified  as  the  layer  depth,  there 
is  a  natural  hesitation  to  proceed  with  sound  field 
calculations  based  on  the  nMinear  profile  mouel. 

Other  choices  of  velocity  profile  models,  complete 
with  solutions  to  the  reduced  wave  equation  that  are 
expressible  in  terms  of  well-known  tabulated  functions, 
are  available.  For  example,  the  reader  is  referred  to 
Murphy  and  Davis  (1974),  wherein  they  discuss  the  use 
of  Weber  functions  when  n2(z)  is  nearly  parabolic. 
The  immediate  question  is,  just  how  many  such 
possible  approaches  are  needed  to  address  the  set  of 
all  possible  profile  configurations? 

There  is  strong  incentive  to  give  serious  considera¬ 
tion  to  propagation  models  that  provide  accurate 
solutions  for  arbitrary  velocity  profile  representations. 
This  point  is  discussed  briefly  as  a  recommendation 
in  Section  V. 

3.  Departure  from  Linearity 

The  illustrations  in  Figure  2  indicate  the  depth  varia¬ 
tion  in  C(Z),  as  well  as  in  n2(Z).  The  variation  in  C 
is  not  quite  linear  in  Z.  A  question  of  interest  is,  what 
impact  does  this  slight  departure  from  linearity 
(if  linearity  is  indeed  the  proper  criterion)  have  on 
various  intermediate  calculations  that  are  involved  in 
producing  an  estimate  of  the  sound  field?  Since  velocity 
profiles  converted  from  in  situ  bathythermographs 
(BTs)  do  not  necessarily  exhibit  a  linear  d'-pth 
dependence,  this  question  should  be  asked  in  reference 
to  the  shape  factor  associated  with  a  given  set  of 
velocity  profile  data.  As  indicated  in  Appendix  F,  the 
n2-linear  model  appears  to  be  a  good  approximation. 

4.  Shape  Factor 

The  previous  discussion  indi.  that  the  curvature 
introduced  by  the  /^-linear  pro-  .  model  is  negligible. 
Suppose  that  actual  profile  data  indicate  that  C5  (Z) 
is  linear  in  depth.  Restricting  attention  to  depths  down 
to  the  layer  depth,  this  profile  model  takes  the  form 

C\Z)  =  C^(l  +  bZ ),  (14) 


Figure  4.  The  n--tinear  and  Ci-linear  velocity  profile  models. 


where  b  =  3 g0/C0  and  g0  is  the  velocity  gradient 
at  the  surface.  Figure  4  illustrates  both  the  nMinear 
and  the  C3-linear  profiles  for  comparison. 

The  CMinear  profile  representation  is  not  as 
popular  as  the  nMinear  model,  even  though  there  are 
surface-duct  environments  for  which  it  would  yield  a 
better  fit.  Unfortunately,  the  corresponding  analytical 
solutions  are  not  expressible  in  simple,  closed  form, 
as  is  the  case  for  the  nMinear  model.  The  corre¬ 
sponding  phase  integral  accedes  to  straightforward 
evaluation,  and  the  resulting  expression  for  the  maxi¬ 
mum  number  of  trapped  modes  is  (see  Appendix  G) 

N  -  V*  +  2A  ( f/ga, )  sin30o/cos20o .  (15) 

By  contrast,  the  corresponding  analysis  for  the 
nMinear  model  leads  to  (see  McGirr,  1983) 

N  =  V*  +  Vi  n-'k0(i»  ZLi/2.  (16) 

For  example,  consider  a  50-m  layer  with  sound 
velocities  of  1500  m/s  at  the  surface  and  1500.9  m/s 
at  the  bottom  of  the  duct.  For  a  constant-gradient 
profile  model  (i.e. ,  C-linear),  these  values  yield  a 
velocity  gradient  of  0.018  sec-1.  For  the  nMinear 
model,  the  parameter  ji  is  2.4  x  10-5  nr1  and  the 
velocity  gradient  as  a  function  of  depth  is  given 
by  dC/dZ  =  Vi{i C3/C2 .  For  the  C3-linear  model, 
the  parameter  b  is  3.6  x  10- 3  nr 1  and  the  velocity 
gradient  as  a  function  of  depth  is  given  by 
dC/dZ  =  VibC\/C2.  The  gradients  evaluated  at  the 
surface  and  at  the  layer  depth  are  summarized  in 
tabular  form  as  follows: 

n2-linear  CMinear 

at  surface  0.0179838  0.0180108 

at  layer  depth  0.0180324  0.0179784 

For  a  frequency  of  5000  Hz,  the  number  of  trapped 
modes  is  7.942  and  7.948  for  the  CMinear  and 
nMinear  models,  respectively.  The  fact  that  the  two 
profile  models  produce  the  same  number  of  modes  does 
not  necessarily  imply  that  the  resulting  sound  field 


estimates  would  be  the  same.  This  agreement,  along 
with  the  good  agreement  between  nMinear  and 
CMinear  expressions  for  the  modal  cycle  range 
(see  Appendix  G),  suggests  that  the  nMinear  profile 
model  is  robust  with  regard  to  the  shape-factor  issue. 

IV.  Coupling  of  Surface  Scattered  Energy 

The  primary  dissipative  contributions  to  total 
surface-duty  propagation  loss  are  (1)  volume  absorp¬ 
tion,  (2)  surface  reflection  loss,  and  (3)  diffractive 
leakage  into  the  negative-gradient  region  below  the 
duct.  Each  of  these  losses  can  be  added  (in  decibel 
space)  to  the  normal-mode  solution  or  they  can  be 
lumped  together  into  a  mode-dependent  attenuation 
coefficient  and  included  in  the  mode  sum.  The  first 
two  forms  of  loss  increase  with  increasing  frequency, 
whereas  the  third  decrease* .  What  is  missing  from  this 
simple  formulation  is  any  prescription  for  interactive 
effects.  A  proper  accounting  of  interaction  between 
surface  scattering  and  leakage  entails  coupling  the 
surface-scattered  energy  directly  into  the  normal-mode 
sum.  Leakage  is  not  the  most  apt  description  of  the 
process,  since  much  of  the  surface-scattered  energy  that 
gets  into  the  below-layer  region  is  by  coupling  between 
trapped  and  untrapped  modes,  as  well  as  by  leaky 
modes. 

When  surface-scattered  energy  is  coupled  into  the 
mode  sum,  the  resulting  sound  field  can  differ 
significantly  from  results  generated  by  a  model  that 
includes  only  a  surface-reflection-loss  term.  That  is, 
when  the  scattering  process  is  treated  as  a  specular 
reflection  process,  there  is  no  mechanism  to  simulate 
the  propagaiion  of  surface-scattered  energy  into  the 
thermocline  region.  For  calm  seas  a  surface-reflection- 
loss  model  can  adequately  account  for  energy  lost  to 
the  scattering  process.  For  fully  developed  seas, 
however,  a  considerable  fraction  of  the  energy  is 
scattered  in  nonspecular  directions. 

These  notions  can  be  put  into  clearer  perspective  by 
giving  some  consideration  to  the  coherent  (/c)  and 
incoherent  or  random  (/,)  components  of  signal 
intensity,  that  is,  /  =  Ic  +  Ir.  The  fraction  of  energy 
in  the  coherent  component,  say  f,  is  given  by 
f  =  /c/(/c  +  Ir).  Let  fs  denote  the  fraction  of  energy 
remaining  in  the  coherent  component  after  reflection 
from  the  surface  and  let  fv  denote  the  fraction 
remaining  after  dissipation  due  to  volume  absorption 
along  the  path  of  propagation.  Then  the  fraction  of 
energy  remaining  in  the  coherent  component  arriving 
at  the  receiver,  say  fR,  is  given  by  fR  =  f?fv,  where 
N  is  the  number  of  surface  reflections. 

Signals  that  propagate  via  a  surface  duct  are 
subjected  to  numerous  encounters  with  the  surface 
boundary.  If  the  surface  is  perfectly  smooth,  i.e., 
fs  =  1,  then  fR  is  essentially  just  fv.  If  the  surface  is 
rough,  however,  then  fs  <  1,  and  for  sufficient 
separation  between  source  and  receiver,  fR  0.  This 


condition  indicates  that  a  substantial  fraction  of  energy 
propagates  incoherently,  or  in  terms  of  the  scattering 
process,  a  significant  amount  of  energy  scatters  off  the 
surface  in  nonspecular  directions.  Thus,  energy  carried 
by  low-order  modes  (small  grazing  angles  for  ray 
equivalents)  can  scatter  off  a  rough  surface  with  steeper 
angles.  If  the  angles  are  steep  enough,  energy  that  is 
ordinarily  trapped  in  the  duct  escapes  into  the  below- 
layer  region.  Most  mode  codes  do  not  account  for 
incoherently  scattered  energy. 

Under  the  proper  conditions,  then,  not  only  can  the 
signal  coherency  (fR/[  1  -  fR])  of  the  trapped  modes 
degrade,  but  a  significant  proportion  of  the 
incoherently  scattered  energy  can  leak  out  of  the  duct. 
For  example,  a  shallow  duct  (small  barrier)  where  both 
in-layer  and  below-Iayer  gradients  are  weak  (small 
interface  reflections)  and  a  rough  sea  surface 
(incoherent  scattering)  provide  conditions  conducive 
to  enhanced  cross-layer  detection  performance  at  the 
expense  of  in-layer  performance.  This  circumstance 
conforms  with  the  situation  described  in  the  introduc¬ 
tory  remarks. 

A  consideration  of  some  importance  is  that  bound¬ 
ary  scattering  is  not  an  inherent  feature  of  wave 
solutions.  This  point  is  often  missed  by  casual  users 
of  propagation  .odes.  What  must  be  kept  in  mind  is 


that  boundary  reflection  loss  models  account  for  the 
fraction  of  energy  that  scatters  only  in  the  specular 
direction.  Bucker  (1980)  takes  surface-scattered  energy 
into  account  by  adding  randomly  phased  scattering 
terms  to  the  coherent  mode  sum.  He  essentially  uses  • 
ray  theory  in  combination  with  mode  theory  to  couple 
surface-scattered  energy  into  the  total  sound  field 
solution.  The  additional  mode-to-ray  and  ray-to-mode 
interaction  terms  are  incoherently  summed  to  obtain 
the  total  intensity  due  to  the  scattered  field.  The 
scattering  (or  interaction)  integrals  developed  by  • 
Bucker  are  presented  in  Appendix  H. 

Bucker  (1980)  demonstrates  the  theory  summarized 
in  Appendix  H  by  comparing  experimental  data  to 
calculations  made  both  with  and  without  the  scattering 
integrals.  The  propagation  loss  calculations  were  made 
using  a  Pedersen-Gordon-type  two-layer,  normal-mode  • 
model.  The  model  inputs  along  with  the  split-beam  and 
surface-reflected  rays  are  presented  in  Figure  5. 

Figure  6  illustrates  two  propagation  loss  vs.  range 
curves.  The  solid  and  dashed  curves  correspond, 
respectively,  to  calculations  made  both  with  and  with¬ 
out  the  scattering  integrals.  The  differences  between  • 
the  two  curves  do  not  become  significant  until  the 
source-receiver  separation  range  reaches  about  7  kyd. 

Beyond  this  range  the  peaks  in  the  multipath  beat 


Figure  6.  Propagation  loss  curves  for  cross-duct  case.  Solid  and  dotted  curves  correspond,  respec¬ 
tively,  to  calculations  made  both  with  and  without  Bucker’s  scattering  integrals. 
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structure  are  1-2  dB  more  optimistic  for  the  solid  curve. 
The  interference  nulls  for  this  curve  are  noticeably  less 
extreme.  This  latter  feature  is  a  characteristic  directly 
attributable  to  the  contribution  of  incoherently  scat¬ 
tered  energy. 

Similar  comparisons  for  the  in-layer  case  are  not 
available,  although  Bucker  (1980)  notes  that  the  differ¬ 
ences  are  not  as  dramatic  for  this  case.  He  further  notes 
that  the  major  differences  occur  in  the  interference 
nulls  where  incoherently  scattered  energy  tends  to  fill 
in  the  sound  field. 

The  interpretation  of  these  results  with  regard  to  the 
,  main  issue— anomalous  SHARPS-generated  surface- 
duct  predictions— is  that  even  though  calculations 
based  on  the  scattering  integrals  yield  more  realistic 
cross-layer  predictions,  the  contribution  of  incoherently 
scattered  energy  is  not  significant  enough  to  explain 
the  anomalous  cross-layer  detection  ranges,  as  reported 
in  the  message  traffic  received  by  FNOC.  Indeed,  the 
differences  between  in-layer  and  cross-layer  predictions 
reported  by  Bucker  (1980)  are  on  the  order  of  2  or 
3  dB — as  determined  from  the  interference  peaks — 
for  ranges  beyond  about  10  kyd.  Thus,  there  appears 
to  be  no  basis  to  support  the  anomalous  cross-layer 
detection  ranges  as  purportedly  predicted  by  SHARPS. 

V.  Computational  Artifacts 

A  problem  that  is  difficult  to  overcome  when  using 
hybrid  solutions  is  the  appearance  of  discontinuities 
that  typically  occur  as  the  method  of  solution  switches 
from  one  method  to  another.  If  the  decision  tree  that 
determines  the  method  of  solution  is  properly  designed 
and  correctly  implemented,  the  magnitude  of  the  jump 
discontinuities  should  be  almost  negligible.  In  Active 
RAYMODE  the  switch  from  the  normal-mode  sum¬ 
mation  method  to  one  of  the  special  integration 
methods  occurs  when  the  number  of  modes  for  a  given 
k-segment  (or  duct)  exceeds  10. 

To  get  an  idea  of  how  serious  this  problem  is  with 
the  Active  RAYMODE  model,  consider  the  following 
example.  The  velocity  profile  is  illustrated  in  Figure  7. 
The  main  interest  here  is  in  the  surface-duct  portion 
of  the  propagation  channel.  Hence,  the  lower  part  of 
the  profile  consists  of  negative  gradients  from  the  layer 
depth  down  to  a  high-loss  bottom  at  a  depth  of  1000  m. 
Both  source  and  receiver  are  in  the  surface  duct  at 
depths  of  6  m  and  18  m,  respectively  .  A  zero  wind 
speed  was  selected  to  discourage  contributions  from  the 
surface  loss  model.  Similarly,  a  high-loss  bottom  was 
specified  to  diminish  contributions  from  bottom- 
reflected  paths. 

Active  RAYMODE  was  executed  at  several  frequen¬ 
cies  in  the  neighborhood  of  2300  Hz  to  expose  any 
discontinuity  that  might  arise  as  a  consequence  of  using 
distinct  methods  of  calculation.  At  this  frequency  the 


number  of  modes  should  just  exceed  10.  That  is,  with 
C0  =  1525  m/s,  ga  =  0.018  s'\  and  ZL  =  100  m, 
the  frequency  corresponding  to  N  modes  is  given  by 
(see  Appendix  G) 

/  =  235.4  (N  -  Vz),  (17) 

which  for  N  =  \0  gives  f  =  2295.2  Hz.  Thus,  there 
should  be  a  discontinuity  at  some  frequency  in  the 
vicinity  of  2300  Hz.  The  scheme  used  in  RAYMODE 
to  determine  the  number  of  modes  trapped  in  a  duct 
is  somewhat  unorthodox  (as  discussed  in  Davis  and 
Council,  1985b),  so  isolating  the  particular  frequency 
at  which  the  switch  occurs  from  the  mode-sum  method 
to  one  of  the  integration  methods  requires  a  trial-and- 
error  search.  A  few  runs  at  several  frequencies  in  the 
vicinity  of  2300  Hz  revealed  that  the  expected  jump 
actually  occurs  at  a  frequency  between  2249  Hz  and 
2250  Hz.  The  discontinuity  is  illustrated  in  Figure  8, 
where  incoherent  transmission  loss  (TL)  is  plotted 
against  frequency  at  a  range  of  5  kyd.  The  data 
shown  are  for  frequencies  varying  from  2200  Hz  to 
2300  Hz  in  10-Hz  steps,  except  for  frequencies  between 
2240  Hz  and  2260  Hz  where  the  step  size  has  been  cut 
down  to  1  Hz.  There  is  a  jump  discontinuity  of  about 
2  dB  between  2249  Hz  and  2250  Hz. 

The  impact  of  switching  from  one  computational 
method  to  another  is  just  as  evident  when  transmis¬ 
sion  loss  vs.  range  curves  are  examined  over  a  band 
of  frequencies  on  either  side  of  the  switching  frequency. 
Figure  9  presents  a  plot  of  two  incoherently  summed 
TL  vs.  range  curves*  for  frequencies  within  a  50-Hz 
neighborhood  on  either  side  of  the  break  frequency. 


•Coherently  summed  TL  is  not  offered  as  an  option  in  Active 
RAYMODE. 
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Figure  8.  Discontinuity  caused  by  switching  from  the  mode  summation  method  to  the  original 
RA  YMODE  integration  method. 


Figure  9.  Transmission  loss  vs.  range  curves  generated  by  Active  RAYMODE  for  the 
frequency  band  less  than  2249  Hz  (symbols)  and  for  the  frequency  band  greater  than  2250 
Hz  (solid). 


According  to  the  documented  decision  tree  (sum¬ 
marized  in  Section  II),  another  discontinuity  is  likely 
to  occur  for  some  frequency  between  2250  Hz  and 
3000  Hz  where  there  is  another  switch  from  the  original 
RAYMODE  integration  (ORI)  method  to  the  special 
HFI  method.  The  only  circumstance  for  which  this 
switch  does  not  take  place  is  when  the  method  of  solu¬ 
tion  switches  directly  from  the  mode-sum  method  to 
the  HFI  method.  A  brief  search  for  the  ORI  ■*  HFI 
switch  revealed  that  the  switch  occurs  at  2752  Hz.  This 
discontinuity  is  illustrated  in  Figure  10.  Here  again  the 
jump  discontinuity  is  on  the  order  of  2  dB. 

The  discontinuity  due  to  switching  from  the  mode- 
sum  method  to  the  ORI  method  also  yields  a  sudden 
jump  for  slight  changes  in  layer  depth.  That  is,  zL, 
corresponding  to  N  modes,  is  given  by 


Substituting  the  values  given  for  c0  and  g0  and  setting 
f  =  2249  Hz  gives,  approximately, 


(19) 


which  for  N  —  10  yields  zL  —  100  m.  Thus,  there 
should  be  a  discontinuity  observed  for  some  layer  depth 
in  the  vicinity  of  100  m.  An  examination  of  this 
particular  type  of  discontinuity  is  not  pursued  here, 
since  it  is  of  the  same  origin  as  the  frequency  discon¬ 
tinuity  illustrated  in  Figure  8. 

The  size  of  jump  discontinuities  of  this  sort  may  not 
be  of  devastating  significance  for  many  cases  of 
interest.  A  cause  of  concern,  however,  is  that  occasions 
can  arise  when  a  mere  2  dB  can  make  the  difference 
between  detection  and  nondetection.  Note  that  the 
slopes  in  the  curves  of  Figure  9  are  not  particularly 
steep  for  ranges  greater  than  about  5  kyd.  A  simple 
(one-way)  figure-of-merit  (FOM)  analysis  reveals  that 
for  an  80-dB  FOM  the  low-frequency  curve  indicates 
a  detection  range  of  about  25  kyd,  whereas  the  high- 
frequency  curve  indicates  that  detections  can  be 
expected  to  cease  for  ranges  greater  than  about  18  kyd. 
There  are  probably  scenarios  (both  active  and  passive) 
for  which  a  discrepancy  of  7  kyd  can  make  a  difference. 


VI.  Summary  Remarks 

The  surface-duct  propagation  loss  prediction 
capability  of  the  Active  and  Passive  RAYMODE 
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Figure  10.  Discontinuity  caused  by  switching  from  the  OR I  method  to  the  HFl  method. 


models  has  beep  qualitatively  assessed;  This  review  was 
motivated  by  reports  (Naval  Messages)  of  questionable 
in-layer  and  cross-layer  target  detection  ranges 
predicted  by  SHARPS  under  certain  well-ducted  prop¬ 
agation  conditions.  Reports  of  this  sort  are  expected 
no  matter  what  model  is  used.  The  largest  variations 
in  sound-velocity  gradients  are  observed  in  the  upper 
portion  of  the  velocity  profile;  therefore,  detection 
range  prediction  errors  are  expected  to  be  greater  for 
near-surface  propagation  paths. 

A.  Conclusions 

The  review  of  RAYMODE  model  physics,  (Section 
II  and  attendant  appendices),  reveals  that  calculations 
pertaining  to  surface-duct  conditions  can  be  performed 
by  one  or  more  methods.  Each  of  the  RAYMODE 
wave-field  integrals  can  take  the  form  of  an  approx¬ 
imate  normal-mode  expansion  or  a  generalized  WKB 
multipath  expansion.  If  the  number  of  propagating 
modes  is  less  than  10,  the  dominant  method  is  an 
approximate  normal-mode  sum.  The  approximations 
that  lead  to  the  RAYMODE  normal-mode  sum  render 
the  approach  essentially  useless  for  surface-duct  cal¬ 
culations.  The  first  approximation  made  is  in  the 
expression  for  the  Wronskian,  where  standard  WKB 
forms  are  used  in  place  of  the  generalized  WKB  ampli¬ 
tude  functions.  The  error  introduced  by  this  particular 
approximation,  however,  should  be  negligible  for  most 
cases  except,  perhaps,  for  surface  ducts.  The  next 
approximation  made  is  in  assuming  real  eigenvalues. 
To  properly  apply  mode  theory  to  surface  ducts,  the 
inclusion  of  modal  attenuation,  due  mainly  to  the  ratio 
of  in-layer  to  below-layer  velocity  gradients,  is  crucial. 

When  the  RAYMODE  wave-field  integrals  are 
evaluated  using  generalized  WKB  expansions,  there 
are  three  distinct  algorithms  from  which  one  is  selected 
according  to  a  decision  tree  (see  Section  II).  Each  of 
these  algorithms  solves  an  integral  of  Fourier  type.  Two 
of  these  methods,  the  ORI  method,  and  the  HFI 
method,  are  reviewed  herein.  Although  stationary- 
phase  (saddle-point)  procedures  are  applied  in  both  of 
these  methods,  in  neither  case  is  the  stationary-phase 


formula  relied  upon.  Instead,  the  integrals  are 
evaluated  over  a  restricted  domain  using  asymptotic 
expansions.  Moreover,  the  HFI  method  includes  a 
special  algorithm  for  the  evaluation  of  a  field  point 
in  the  vicinity  of  a  caustic.  The  decision  rule  encom¬ 
passes  the  approximate  normal-mode  sum  as  well;  thus, 
the  RAYMODE  decision  process  effectively  partitions 
the  wave-number  integration  into  low-  and  high- 
frequency  components.  Since  this  partitioning  is 
accomplished  for  a  given  duct,  there  is  always  the 
possibility  that  several  of  the  RAYMODE  methods  are 
exercised  during  a  single  model  execution. 

The  magnitude  of  error  that  might  be  introduced  by 
the  ORI  and  the  HFI  methods  is  difficult  to  assess  in 
general.  The  accuracy  of  these  methods  is  closely 
associated  with  the  validity  of  WKB  solutions  to  the 
depth-separated  wave  equation.  The  local  sound  field 
is  well  approximated  by  WKB  solutions  for  certain 
depth  intervals  referred  to  as  “allowed”  regions.  These 
methods  may  be  inaccurate  in  complementary  intervals, 
referred  to  as  “forbidden”  regions.  In  an  attempt  to 
improve  standard  WKB  solutions,  the  model  developer 
uses  a  generalized  WKB  solution,  which  tends  to  degen¬ 
erate  to  the  standard  WKB  form  only  when  the  field 
point  is  far  removed  from  a  turning  point.  The 
generalized  WKB  amplitude  functions  are  determined 
iteratively.  The  only  other  instance  when  the  standard 
WKB  form  is  relied  upon,  is  in  arriving  at  an 
approximate  expression  for  the  Wronskian.  This 
approximation  impacts  all  of  the  RAYMODE  methods 
reviewed  herein.  A  previous  RAYMODE  evaluation 
conducted  by  Bartberger  (1981)  indicates  that  the  error 
introduced  by  the  special  integration  methods  has  its 
major  impact  on  coherently  summed  results,  being 
somewhat  mitigated  when  the  results  are  incoherently 
summed.  The  general  requirements  justifying  a 
stationary-phase  analysis  are  probably  not  met  at  low 
frequencies,  i.e.,  phase  oscillations  may  not  be  rapid 
enough  to  effect  the  required  cancellations.  As  far  as 
the  “extension”  of  this  method  to  include  those  cases 
when  the  point  of  stationary  phase  falls  outside  the 
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limits  of  integration  is  concerned,  the  procedure  used 
in  RAYMODE  appears  to  have  been  chosen  out  of 
convenience,  i.e.,  so  that  the  Fresnel  integral  algorithms 
can  be  used.  The  specific  steps  taken  by  the  model 
developer  in  applying  stationary-phase  methods  are 
not  universally  accepted  as  being  the  most  effective 
(see  Weinberg,  1981). 

The  surface  reflection  loss  submodel  is  reviewed  in 
Section  II.  An  estimate  of  surface  loss  (SL)  is  deter¬ 
mined  from  the  sum  of  two  independent  terms,  SI, 
and  SL2.  The  evaluation  by  Keenan  (1983)  of  the 
SUBCOM  surface  loss  model  (the  one  used  in 
RAYMODE)  shows  that  the  SI,  term  represents 
incoherently  scattered  energy  characteristic  of  a 
Kirchhoff  solution.  Keenan  further  notes  that  this  term 
is  dominant  for  frequencies  below  about  1  kHz.  Since 
incoherent  scattering  should  dominate  the  solution  only 
for  large  values  of  the  Rayleigh  parameter,  Keenan 
concludes  that  the  SI,  term  is  not  intuitively  satisfac¬ 
tory  as  an  estimator  of  rough  surface  loss. 

The  SI 2  term  has  no  angle  dependence,  and  it  is 
probably  supposed  to  represent  the  coherent  compo¬ 
nent  of  energy  scattered  at  an  average  angle 
characteristic  of  surface-duct  and  convergence  zone 
paths.  Speculation  of  this  sort  is  interesting,  although 
not  particularly  productive.  Until  the  model  developer 
subjects  his  derivation  to  peer  review,  however,  any 
critique  of  the  model  physics  amounts  to  little  more 
than  speculation.  Results  generated  by  this  model  were 
evaluated,  along  with  results  generated  by  several  other 
surface  loss  models.  Comparisons  are  included  in  a 
report  that  presents  the  findings  and  recommendations 
of  a  surface  loss  model  working  group  (Eller,  1984). 
Interestingly,  the  surface  loss  model  recommended  by 
this  working  group  has  not  been  adopted. 

A  special  surface-duct  leakage  routine,  included  in 
a  previous  version  of  Passive  RAYMODE,  is  con¬ 
spicuously  missing  from  the  present  version.  A  review 
of  the  corresponding  computer  code  reveals  at  least 
one  error,  along  with  expressions  that  bear  a  curious 
resemblance  to  those  normally  associated  with  WKB 
phase  integral  methods.  Corrections  could  be  made 
quite  readily,  although  suggestions  offering  more 
promise  are  made  in  the  recommendations  section. 

Before  a  two-layer  surface-duct  submodel  is  decided 
upon,  serious  consideration  should  be  given  to  the  con¬ 
sequences:  namely,  its  inflexibility  with  regard,  not  only 
to  velocity  profile  fit,  but  to  its  inability  to  address 
more  general  near-surface  profile  configurations.  If 
such  a  submodel  is  selected,  there  are  ways  to  improve 
computational  efficiency  (see  Section  III). 

The  omission  of  any  means  to  couple  surface- 
scattered  energy  into  the  sound  field  is  not  a  deficiency 
peculiar  to  the  RAYMODE  codes.  Indeed,  few  models 
address  this  important  feature.  For  the  example 
discussed  in  Section  IV,  the  inclusion  of  Bucker’s  scat¬ 
tering  integrals  can  bring  the  interference  peaks  up 


(reduce  the  loss)  by  2  dB  for  the  longer  ranges.  Just 
as  importantly,  the  interference  nulls  tend  to  get  washed 
out.  For  detection  range  calculations  based  on  coherent 
signal-to-noise  ratios,  the  contribution  of  incoherently 
scattered  energy  would  be  significant,  but  only  if  the 
correct  statistical  distributions  are  considered  (as 
discussed  in  Pedersen  and  McGirr,  1983).  If  only 
averaged  signal-to-noise  ratios  are  considered,  the 
significance  is  less  dramatic.  Consideration 
should  be  given  to  including  the  scattering  integrals, 
since  their  omission  can  yield  slightly  pessimistic  cross- 
layer  predictions. 

Discontinuities  that  occur  when  the  method  of  solu¬ 
tion  switches  from  one  form  to  another  appear  to  be 
on  the  order  of  a  couple  of  decibels,  at  least  as  deter¬ 
mined  by  the  brief  assessment  presented  herein. 
Numerical  artifacts  of  this  sort  fall  into  the  category 
of  nuisances  and  are  typical  of  hybrid  methods.  To 
the  extent  that  the  size  of  these  discontinuities  does  not 
exceed  about  2  dB,  they  are  not  considered  to  be  major 
deficiencies. 

In  any  case,  none  of  the  Active  RAYMODE  model 
deficiencies  discussed  herein  is  likely  to  cause 
anomalous  cross-layer  detection  ranges,  such  as  those 
purportedly  predicted  by  SHARPS.  The  very 
appearance  of  anomalous  surface-duct  detection  ranges 
generated  by  SHARPS  naturally  raises  the  question, 
what  set  of  surface-duct  conditions  could  possibly 
produce  greater  detection  ranges  for  cross-layer 
geometries  than  for  in-layer  geometries?  For  a  well- 
defined  classic  surface  duct,  an  answer  is  not 
immediately  forthcoming.  If  care  is  not  taken  in 
identifying — acoustically— what  does  or  does  not 
constitute  a  well-defined  surface  duct,  then  interpreta¬ 
tions  of  the  resulting  acoustic  predictions  are 
more-or-less  at  the  mercy  of  possibly  invalid  interpreta¬ 
tions  of  the  controlling  velocity  profile.  This  latter 
possibility  essentially  reflects  the  essence  of  the 
so-called  sonic-layer-depth  (SLD)  problem. 

As  a  case  in  point,  consider  the  following  example. 
Suppose  the  input  velocity  profile  (as  generated  by  a 
pre-SHARPS  program)  takes  the  form  illustrated  in 
Figure  11.  The  near-surface  portion  of  this  profile 
consists  of  an  isovelocity  duct  overlaying  a  weak 
gradient  layer.  Even  though  the  SLD  for  this  case 
would  be  reported  in  the  SHARPS  message  as  50  m, 
the  first  two  layers  actually  comprise  the  in-layer 
portion  of  a  weak  surface  duct.  The  difference  in  the 
velocity  gradients  of  these  two  layers  is  acoustically 
negligible.  The  major  gradient  discontinuity  occurs 
between  the  second  and  third  layers;  therefore,  the 
bottom  of  the  in-layer  portion  of  the  near-surface 
channel  actually  corresponds  to  150  m. 

For  demonstration,  however,  let  the  upper  two  layers 
be  arbitrarily  identified  as  the  in-layer  and  below-layer 
components  of  the  surface  duct,  which  would  be  the 
logical  interpretation  formulated  by  a  recipient  of  a 
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SONAR  DEPTH  =  6  m 

TARGET  DEPTH  =  18  m/100  m 

FREQUENCY  =  3000  Hz 
WIND  SPEED  =  Okt 

1535  m/s 


Figure  11.  Sound  velocity  profile  for  an  SLD  test  case. 


SHARPS  message  claiming  a  50-m  SLD.  Let  the  sonar 
and  target  depths  be  6  m  and  18  m,  respectively,  for  the 
“in-layer”  case,  and  6  m  and  100  m,  respectively,  for 
the  “cross-layer”  case.  Figure  12  presents  (a)  in-layer 
and  (b)  cross-layer  propagation  loss  vs.  range  curves 
for  a  frequency  of  3000  Hz.  Not  surprisingly,  the 
cross-layer  curve  is  more  optimistic  than  the  in-layer 
curve.  By  virtue  of  their  weak  gradients,  the  first  two 
layers  of  this  profile  do  not  constitute  a  well-defined 
surface  duct,  but  because  of  the  way  SLD  is  defined, 
curves  (a)  and  (b)  of  Figure  12  would  be  interpreted 
as  predictions  for  in-layer  and  cross-layer  geometries. 

B.  Recommendations 

The  most  important  deficiency  in  the  RAYMODE 
codes  that  needs  to  be  rectified  is  the  omission  of  any 
means  to  account  for  interduct  coupling.  This  problem, 
in  general,  is  not  easily  corrected  without  recourse  to 
rigorous  normal-mode  theory;  even  then,  numerical 
problems  can  arise  during  attempts  to  And  eigenvalues 
in  interduct  situations.  The  easiest  solution  is  to  add 
a  special  surface-duct  module,  and  simply  ignore  other 
interduct  coupling  problems.  One  possibility  is  to  resur¬ 
rect  the  surface-duct  submodel  from  a  previous  ver¬ 
sion  of  Passive  RAYMODE.  A  brief  discussion  of  this 
submodel  is  presented  in  the  document  of  Davis  and 
Council  (1985b),  as  well  as  in  Section  II.  Some 
additional  work  on  this  submodel  is  indicated. 
ASTRAL*  already  has  such  a  module  (Ryan,  1980; 
White,  1988);  therefore,  it  could  be  incorporated  rather 
easily  into  the  RAYMODE  codes.  Inasmuch  as 


•ASTRAL:  ASEPS*  Transmission  Loss  (Spofford,  1979) 
'ASEPS:  Automated  Signal  Excess  Prediction  System 
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ASTRAL  is  one  of  the  standard  range-dependent 
propagation  loss  models,  this  approach  would 
automatically  satisfy  the  OPNAV  commonality 
requirement. 

The  ASTRAL  surface-duct  module  does  not  deal 
with  energy  that  is  incoherently  scattered  by  the 
surface.  To  properly  address  surface-scattered  energy, 
this  module  would  have  to  be  upgraded.  The  incor¬ 
poration  of  Bucker’s  (1980)  scattering  integrals  is  easy 
conceptually,  but  could  be  costly.  Actually  there  are 
three  costs  to  consider:  (1)  the  cost  of  implementation, 
and  once  implemented,  the  costs  in  (2)  storage  and  in 
(3)  run-time  that  would  be  incurred  operationally. 

Instead  of  using  a  two-layer  surface-duct  model, 
serious  consideration  should  be  given  to  using  a 
normal-mode  model  that  can  handle  up  to,  say,  five 
layers.  In  this  way,  accurate  predictions  can  be  made 
for  all  important  near-surface  velocity  profile  con¬ 
figurations.  The  major  objections  to  this  approach 
naturally  pertain  to  excessive  run-time  and  storage 
requirements,  but  the  rate  of  advances  presently  being 
made  in  computer  technology  suggests  that  these 
requirements  are  achievable  in  the  very  near  future  for 
both  ashore  and  afloat  systems.  Moreover,  as  available 
computing  resources  become  more  powerful,  the 
emphasis  in  modeling  upgrades  should  focus  on 
improving  the  quality  of  predictions  vice  quantity. 

This  recommendation  is  not  entirely  new  (McGirr, 
1983)  although  a  full-wave  treatment  is  recommended 
instead  of  an  approximate  approach  based  on  WKB 
phase-integral  methods.  Consider  the  expression  that 
gives  the  (far-field)  acoustic  pressure  for  a  receiver  at 
depth  z  and  range  r,  i.e., 

p(r,z)  •  C  I,  uH(zJ  u„(z)  expiikHr)/(knr)'/1,  (20) 
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where  z0  is  the  source  depth,  and  the  un  are  normal¬ 
ized  modal  depth  functions.  The  kn  (eigenvalues)  are 
functions  of  frequency  and  the  environmental-acoustic 
input  profile  (velocity  profile,  bottom  depth,  bottom 
loss,  and  surface  loss),  but  they  are  independent  of 
source  depth  and  receiver  depth.  The  u„(z)  are  func¬ 
tions  of  the  kn.  Thus,  there  is  a  natural  hierarchical 
ordering  of  the  required  computations  for  any  given 
run. 

Suppose  that  predictions  are  required  for  a  single 
source  depth  and  three  receiver  depths  at  some 
frequency,  say,  f0,  for  which  there  are  N  modes.  First 
the  N  eigenvalues  are  determined  and  stored.  Then, 
the  depth  functions  {«„(*<,)}  are  computed  and 
stored.  Then,  three  sets  of  /V depth  functions  {ujiz)}, 
j  =  1,2  and  3,  are  computed  and  stored.  For  each 
receiver  depth  (z  ),  a  transmission  loss  vs.  range  curve 
can  be  determined  from 

TL  10  log  |pp*  .  (21) 

where,  for  long  ranges,  the  normalized  acoustic 
pressure  is  given  by 

p  =  2  An  u„(z)  exp (iknr)/y/r,  (22) 

and  where  An  =  C  u„(zj/  kn.  The  expression  for  p 
has  been  put  into  this  form  to  emphasize  that  the  An 
are  simply  retrieved  from  storage  during  this  stage  of 
computation.  The  number  of  terms  in  the  sum  is 
actually  a  function  of  range.  That  is,  modal  attenua¬ 
tion  increases  with  range,  therefore,  fewer  terms  are 
needed  as  range  increases.  A  relatively  conservative 
criterion  used  in  determining  when  to  stop  the  sum¬ 
mation  entails  comparing  the  last  four  terms  to  the 
cumulative  sum.  If  the  relative  contribution  is  less  than 
0.0001,  then  the  summation  is  stopped. 

The  computational  strategy  described  above  can  be 
exploited  even  more  aggressively  in  the  design  of 
prediction  systems  dedicated  to  the  generation  of  data 
in  support  of  predeployment  acoustic  assessments. 
Assessments  of  this  sort  are  typically  based  on  archived 
data.  Standard  Navy  data  bases  could  be  expanded  to 
include  precalculated  intermediate  modal  quantities 
tailored  to  specific  system  parameters  and  most- 
probable  target  depths.  This  procedure  is  similar  to 
the  one  discussed  in  Section  II.  To  incorporate  updated 
oceanographic  data,  however,  would  require  adding 
a  perturbation  scheme,  or,  using  a  prediction  system 
such  as  WRAP,  which  is  presently  under  development 
by  Kuperman  et  al.  (1986)  and  Porter  et  al.  (1987). 

Finally,  insofar  as  anomalous  SHARPS-generated, 
cross-layer  detection  ranges  are  concerned,  this  cir¬ 
cumstance  is  more  than  likely  a  consequence  of  the 
long-standing  SLD  problem,  as  mentioned  in  the  intro¬ 
duction  and  discussed  in  the  conclusions.  This  problem 


is  also  d^cussed  by  McGirr  (1983)  and  by  Renner  and 
Kirby  (1983).  Since  the  anomalous  SHARPS  predic¬ 
tions  are  in  contrast  with  those  generated  by  shipboard 
prediction  systems,  (the  latter  results  are  purportedly 
the  more  reasonable)  the  OPNAV  imposed  com¬ 
monality  requirement  needs  to  be  extended  to  include 
any  and  all  peripheral  software  and  data  bases  that 
can  impact  predictions  generated  by  standard  acoustic 
models.  Thus,  as  a  final  recommendation,  the  organi¬ 
zation  responsible  for  configuration  control  of  Fleet 
acoustic  prediction  systems  should  ensure  that  all  such 
software  packages  and  data  bases  essentially  produce 
identical  results.  For  the  case  at  hand,  this  action  entails 
replacing  the  BT-conversion  software  presently  used 
at  FNOC  with  the  appropriate  software  that  is 
used  by  (standard)  shipboard  prediction  systems. 
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Normal-mode  and  Multipath  Expansion  Forms 


k 


k 


k 


k 


k 


i 


> 


i 


The  review  presented  herein  focuses  on  those  methods  that  are  most 
apt  to  have  some  bearing  on  surface  duct  propagation  loss  predictions. 
Since  the  propagation  model  used  in  Active  RAYMODE  is  essentially  a 
subset  of  Passive  RAYMODE,  most  of  the  discussion  presented  in  the 
following  paragraphs  serves  as  a  review  of  those  methods  used  in  Passive 
RAYMODE  that  can  impact  predictions  for  surface  duct  environments. 
Indeed,  at  the  outset  of  this  task,  the  hypothesis  was  formulated  that  the 
special  surface  duct  treatment  found  in  a  previous  version  of  Passive 
RAYMODE  can  be  used,  with  corrections  and  suitable  modifications,  to 
couple  the  surface  scattered  field  into  the  mode  sum,  thereby  yielding 
better  sound  field  estimates  for  both  in-layer  and  cross-layer  geometries. 

Passive  RAYMODE  is  an  interesting  model  because  of  the  several 
methods  that  are  applied  to  the  problem  of  solving  the  wave  equation.  Not 
all  of  these  methods  are  applicable  to  the  surface  duct  problem  of 
interest  here,  and  hence  attention  focuses  on  those  methods  that 
potentially  have  some  bearing  on  the  problem  at  hand.  A  convenient 
starting  point  for  this  review  is  the  wave  equation,  which,  under  the 
assumption  of  azimuthal  symmetry  and  range-invariant  horizontal 
stratification  of  the  water  column,  can  be  expressed  as 


ILfril) i  a2Y  _  j_  a2Y 

r  dr'  dr  y  *  $ z 2  "  c*  IP  ' 


(A  I ) 


where  r  is  range,  z  is  depth,  t  is  travel  time,  c  is  the  velocity  of  sound 
treated  as  a  function  of  depth  only,  and  \y  is  the  velocity  potential. 
Assuming  a  point  harmonic  source  with  exp(icot)  time  dependence,  the 
velocity  potential  can  be  expressed  as 

y  «  eio)t<I >(r,z),  (A2) 

where  i  «  V-1 ,  o>  -  27tf,  and  $  is  a  potential  function  dependent  on  the 
spatial  variables  r  and  z.  Substitution  of  Eq.  (A2)  into  Eq.  (A1)  yields 


r  dr 'dr  )  $z2 


.2 
d  * 


(A3) 
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Since  surface  duct  propagation  is  the  main  focus  of  this  report,  there 
is  no  need  to  include  effects  due  to  variations  in  density,  and  hence  the 
density  can  be  taken  as  constant  (i.e.,  p  =  1).  Thus,  the  acoustic  pressure  p 
is  simply 

p  -  icocD.  (A4) 

Expressing  O(r.z)  as  the  product  of  R(r),  a  function  of  range  only,  and  Z(z ), 
a  function  of  depth  only,  leads  to  the  following  pair  of  equations 


and 


*d*R  rdR 

dr*  dr 

2  2 

+  kr  R  =  0, 

(A5) 

ii  *  (s.2 

iz1  \c2 

-k2)z  =0, 

(A6) 

where  k  is  the  horizontal  wave  number.  Eq.  (A5)  is  recognized  as  a  form  of 
Bessel's  equation  which  has  the  solution 

R(r)  -  J0(kr),  (A7) 

where  J0(kr)  is  the  zero-order  Bessel  function  of  the  first  kind.  The  usual 
practice  in  dealing  with  underwater  sound  problems  is  to  consider  the 
Fourier-Bessel  transform,  that  is 

(1)  °?  (2) 

*  =  J  G(zr,zs,k)H0  (kr)kdk  +  ]G(zr,z  k)H0(kr)kdk,  (A8) 

O  0 

where  J0(kr)  has  been  replaced  by  the  zero-order  Hankel  functions  of  the 
first  and  second  kind  through  the  relation 

J0(kr)  -  i/2  [H00)(kr)  +  H0(2)(kr)].  (A9) 

Since  H^^kr)  ■  -H0<2)(-kr),  Eq.  (A8)  can  be  expressed  in  the  form 

oo 

<D-  J  G(zr,zs,k)  H0(2)(kr)  k  dk 


(A10) 


where  G(zr,zs;k)  is  referred  to  as  the  depth-dependent  Green's  functions 

The  development  from  this  point  follows  that  of  Davis  and  Council 
(1985b),  although  some  liberties  have  been  taken  in  making  slight 
modifications  in  notation.  Let  U(z,k)  and  V(z,k)  be  linearly  independent 
solutions  of  Eq.  (A6),  where  U  satisfies  the  upper  boundary  condition  and  V 
satisfies  the  lower  boundary  condition.  Denoting  source  depth  by  zs  and 
receiver  depth  by  zr,  the  Green's  function  takes  the  form 


G  -  U(zs,k)  V(zr,k)  /  W,  i*  zr  >  zs, 

(All) 

or 

G  «  U(zr,k)  V(zs,k)  /  W,  if  zr  <  zs, 

( A 1  2) 

where  W,  the  Wronskian,  is  given  by 

W  -  U  dV/dz  -  dU/dz  V. 

(A13) 

At  this  stage  the  approach  used  in  P.AYMGDE  is  to  assume  that  U  and  V 
can  be  expressed  in  terms  or  generalized  WKB  solutions  that  take  the 
forms 

U  -  u(z)  exp[  i$(zu,zs)j  +  Ru  u*(z)  exp[-i<|>(zu,zs)]  (A14) 

for  zu  <  zs,  and 

V  .  v(z)  exp[-i4>(z,,zr)]  +  R,  v*(z)  exp[  i<t>(z,,zr)]  (A15) 

for  zr  <  z\.  A  tacit  assumption  made  in  these  expressions  is  that  zs  <  zr, 
although  there  is  no  loss  of  generality  since  their  roles  can  be  reversed  by 
a  simple  interchange  of  subscripts.  Ru  and  R(  are  reflection  coefficients 
that  must  satisfy  certain  conditions  at  the  upper  and  lower  boundaries  or 
turning  points,  respectively.  In  an  attempt  to  improve  upon  standard  WKB 
solutions,  the  amplitude  functions,  u(z)  and  v(z),  are  determined 
iteratively. 


tFor  more  details  see  Sect.  47  of  Brekhovskikh  (1980) 


The  function  appearing  in  the  exponents  is  given  by 


<j>(Zt,Z2) 
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1  Q(z,k)  dz, 
Zl 


where 


Q2(z.k)  -  [<a/c(z)]2  -  k2. 


( A 1  6) 


( A 1  7) 


The  upper  and  lower  turning  points  occur  at  depths  zu  and  Zi  where  the 
function  Q  vanishes. 

The  Wronskian  that  appears  in  the  denominators  of  Eqs.  (All)  and 
(A12)  can  b  •  approximated  by 


where 


W  =  -i2  (1  -  A)  exp{-i0(zu,z,)], 
A  -  RUR,  exp[-i2  4>(zu,z,)]. 


(A18) 

(A19) 


The  approximation  that  leads  to  this  simple  expression  for  W  rests  on  the 
assumption  that  the  WKB  amplitude  functions  can  be  represented  by 
standard  WKB  forms,  i.e.,  u  -  v  2  Q-1/2. 

Using  the  first  term  in  the  asymptotic  expansion  of  H0(2)(kr),  i.e., 


Ho(2)(k0  s  (2/nkr)1/2  exp[-i(kr  -  w/4)], 
the  integral  in  Eq.  (A10)  can  be  expressed  in  the  form 


(A20) 


4  dk 

expH<V„*kr)]  — 


(A21) 
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where  the  Fn  are  given  by 


and 

with  K(kr) 


Fi  -  u  (zs)  v  (zr)  K(kr), 

F2  -  Ru  u*(zs)  v  (zr)  K(kr), 

F3  -  R,  u  (zs)  v*(zr)  K(kr), 

F4  -  RUR,  u*(Zg)  v* (zr)  K(kr), 

(k/2 *r)1/2  and  where  the  \yn  are  given  by 


(A22) 


V,  -  0(ZU,Z|)  -  <D(zu,zs)  -  <>(Zr.Z|), 
v2  -  <t>(Zu ,Z|)  +  4>(ZU,ZS)  -  0(Zr,Z|), 

V3  -  0(zUlZ|)  -  <>(zu,zs)  +  <>(zr,Z|), 

V4  -  <j>(Zy ,Z|)  +  <t> (z (j , zs)  +  <t»(zr,Z|).  (A23, 


Thus,  the  integral,  Eq.  (A10),  has  essentially  been  decomposed  into  four 

wave  field  integrals,  each  of  which  can  be  associated  with  a  ray  path  that 

leaves  the  source  in  the  up  or  down  direction  and  arrives  at  the  receiver  in 
the  up  or  down  direction.  The  particular  identifications  are  given  in  the 
RAYMODE  documents  of  Davis  and  Council  (1985a,b). 

Generally,  U  and  V  assume  different  forms  depending  on  which  side  of 
a  turning  point  depth  the  field  point  depth,  z,  is  located.  If  the  velocity 
minimum  is  at  a  depth  falling  between  the  turning  point  depths,  then  for 
zu  <  z  <  Z|  ("allowed"  region)  these  functions  are  oscillatory,  whereas 
for  z  <  zg  or  z  >  Z|  ("forbidden"  region)  they  take  on  exponentially  decaying 
forms.  In  either  of  these  regions,  the  amplitude  functions,  u(z)  and  v(z), 

reduce  to  the  standard  WKB  form,  i.e.,  ~Q*1/2(z),  only  when  the  field  point 

is  far  removed  from  a  turning  point.  More  generally,  however,  they  are 
determined  iteratively  using  recurrence  relations  (Leibiger,  1971). 


One  of  the  unique  features  of  RAYMODE  is  the  manner  in  which  the 
limits  of  integration  are  handled.  When  the  integration  and  summation 
operations  indicated  in  Eq.  (A21)  are  interchanged,  the  bulk  of  the 
computational  effort  can  be  seen  to  center  on  integrals  of  the  form 


I.  MOO  ej|Wkrl 

-oo 


dk 

1  -  A(k)  ’ 


(A24) 


Leibiger  (1971)  recognized  that  computational  economies  could  be 
realized  by  judiciously  partitioning  the  k  axis,  assumed  to  be  the  real  line 
(•«,«),  into  several  much  smaller  segments.  Thus  the  integration  of  Eq. 
(A24)  is  performed  in  a  piecewise  manner,  that  is, 

oo  kM-«  kn-«  kn+1-<  kM“« 

J  —►  /  =/+/♦•+/  <A25> 

-0°  k„.  .  k_,  . 


The  set  of  k  values  is  first  truncated  to  the  interval  (km,  kM)  which  is  then 
further  subdivided  into  segments  (kn+e,  kn+1_e),  where 

km  ^  kn+e  <  ^n+1-e  ^ 


and  £  is  an  arbitrarily  small  parameter  introduced  to  ensure  the  integrity 
of  the  subsegments.  The  selection  of  km,  kM  and  the  further  segmentation 
of  (km-  I'm)  is  discussed  in  some  detail  in  Davis  and  Council  (1985b),  who 
consider  a  three-duct  profile,  and  in  Deavenport  (1978),  who  considers  a 
classical  deep-ocean  profile  containing  a  surface  duct.  In  his  discussion 
of  integration  limits,  Deavenport  notes  that  (referring  to  a  previous 
version  of  RAYMODE)  surface  duct  situations  are  subjected  to  special 
considerations  based  on  complex  k  values.  These  considerations  are 
reviewed  in  Section  III  (and  Appendix  E). 


For  a  given  k-segment,  say  ka  to  kb,  the  integrals  take  the  form 


_-ll'Kk)+  krl  dk 

I(ka,kb)  -  [  FCk)  9  t  _  A(k) 


(A26) 
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When  the  number  of  k  values  is  not  excessive,*)*  the  integrals  are  evaluated 
using  normal  mode  theory.  The  eigenvalues  required  in  the  mode  sum  are 
determined  by  solving 

1  -  A(km)  -  1  -  RUR,  exp(-i2$m)  -  0,  (A27) 

where  4>m  s  $(zu,Z|),  and  the  turning  point  depths  zu  and  z,  correspond  tc 
mode  m.  Equation  (A27)  can  be  expressed  as 

1  -  IRulelS*u  ll^le1**1  eH2*m=  0.  (A28) 

The  terms  8$u  and  8<J>|  are  phase  changes  that  occur  at  the  upper  and  lower 
turning  points.  Representing  <|>m  as  a  complex  quantity  with  real  part 
denoted  by  and  (small)  imaginary  part  denoted  by  pm,  the  real  and 

imaginary  components  of  Eq.  (A28)  are  given  by 

(real  part)  1  -IRUR,I  e2Pm  cos(  -  2 Tm)  =  o,  (A29) 

and 

(imag  part)  iRyR,!  e2Pm  $1n(  *♦*,  +  *♦,  -2?m)  =  o.  (A30) 

From  the  imaginary  part 

S0U  +  8<>,  -  2  -  2(m-1)«,  m  -  1,  2,  ...  (A31) 

which  implies  that  cos(5$u  +  8$i  -  $m)  ■  1,  so  that,  from  the  real  part, 

IRUR,I  =  B'*fm  .  (A32) 

By  assuming  that  km  -  km  +  iam  and  that  am  (*  lm{km})  contributes 
mainly  to  modal  attenuation,  the  eigenvalues  are  taken  to  be  simple  poles 
of  the  integral  (Eq.  ‘  (A26))  or  the  zeros  of  Eq.  (A27).  Equivalently,  the 
condition  expressed  by  Eq.  (A31)  may  be  used  if  the  mode  index,  m, 


t  the  maximum  number  of  modes  is  currently  limited  to  30  in  Passive 
RAYMODE  and  10  in  Active  RAYMODE 


is  temporarily  treated  as  a  continuous  variable.  That  is,  a  zero  of  1  -  A(k) 
can  be  found  by  solving  d  $m/dm  -  n.  Since 

-  Zl  -1/2 

si*  £mi«V-i£i  dz.  (A33) 

Zg 

and  since  this  integral  is  just  the  range  of  a  ray  proceeding  from  zu  to  Z\, 
corresponding  to  one-half  cycle  for  mode  m,  then 

=  -^[Re(8mV2).  (A34) 

Thus  d  km/drn  -  -2jc/Rc(0m),  where  Rc(0m)  is  the  cycle  range  for  a  ray 
launched  with  angle  8m.  Differentiation  of  Eq.  (A27)  w.r.t.  km  -  taking  Eq. 
(A34)  into  consideration  -  yields 

fj-O-A^-1  «.«„>.  (A35) 

Next  the  phase  factors,  4>m(k)  and  ym(k),  are  approximated  by  low-order 
expansions,  e.g., 


=  4>(  km)  +  ^>m(  ^m)  0*  "  ^m) 

2  <t>(  km)  ’1/2Rc(0m)ctm*  (A36) 

These  approximations  along  with  the  expansion  of  residues  yields  the 
RAYMODE  normal  mode  sum,  that  is 


ReCm) 


-UYCO+  kmr] 

6  8 


(A37) 


where  the  reflection  coefficients  have  been  treated  as  constants. 
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When  the  number  of  modes  is  excessive,  one  of  several  other 
approaches  can  be  used  in  performing  the  integrations  indicated  in 
Eq.  (A26).  Before  one  of  these  other  techniques  is  applied,  however,  Eq. 
(A26)  is  recast  into  "multipath  expansion"  form.  That  is, 


I(kaX)  =  I 

*  b  n«o 


|RyR,|  F  expHfv  +2n$(zu,z.)  ♦  kr]}  dk,  (A38) 


where  [1  -  A(k)]"1  has  been  replaced  by  the  infinite  seriest  I  An(k).  This 
series  is  convergent  for  any  physical  case  of  interest  since  |A(k)|  <  1.  This 
stipulation  on  the  magnitude  of  A  (  »  |RUR ,|)  naturally  excludes  "perfectly 
reflecting"  boundaries,  but,  of  course,  boundaries  of  this  sort  do  not  exist 
in  actual  ocean  environments.  Note  that  Eq.  (A38)  is  but  one  of  four 
integrals  that  comprise  the  total  solution,  and,  in  this  regard,  the  reader 
is  reminded  that  the  forms  taken  on  by  F  and  y  differ  slightly  from  one 
integral  to  the  next.  Using  the  enumeration  scheme  of  Eq.’s  (A22)  and 
(A23),  when  the  reflection  coefficient  magnitudes,  |Ru|n  and  |R(/n,  are 
incorporated  into  the  F-factors,  each  F-factor  has  |RU|  raised  to  the  power 
n  +  p,  where  p«  0,1,0  and  1,  respectively,  and  |R||  raised  to  the  power 
n  +  q,  where  q  •  0,  0,  1  and  1,  respectively. 


t  For  an  elementary  interpretation  of  this  technique  consult  Sect.  35 
in  the  text  by  Brekhovskikh  (1980). 
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Appendix  B 

Special  RAYMODE  Integration  Methods 


The  general  form  of  Eq.  (A38)  is 


l#,V  =  f6p(k)e'1tkr*a<k)1dk.  (01) 

k4 

An  approach  often  used  in  the  evaluation  of  integrals  of  this  form  is  the 
method  of  stationary  phase.  A  point  of  stationary  phase  occurs  at  the 
point  k*,  say,  such  that  the  argument  of  the  exponential  function  -  the 
phase  -  has  zero  slope.  That  is, 

|^[kr  +  a.oo]|  #  =  r+  <x(k#)  =  0,  (B2) 

k 

where  the  prime  indicates  partial  differentiation  with  respect  to  k.  The 
contributions  to  the  integral  from  values  of  k  outside  the  interval  (k*  - 
Ak,  k*  -  Ak)  tend  to  cancel  due  to  oscillations  in  phase,  but  inside  this 
interval  the  value  of  kr  +  a(k)  is  nearly  constant. 

Following  the  development  of  Erd6lyi  (1956),  let  h(k)  -  kr  +  a(k)  so 
that  Eq.  (Bl)  takes  the  form 

kb  -ihfkl 

1=  f  P(k)e  %k.  (03) 

k» 

If  k*  is  a  point  of  stationary  phase,  then 

h(k)  s  h(k*)  +  V2  h"(k*)  (k  -  k*)2.  (B4) 

Now  let  u2  -  h(k)  -  h(k*)  so  that,  by  implication, 

u  -  [h"/2]1/2(k  -  k*),  (B5) 

where  h"  •  h"(k)|k*.  With  this  substitution,  the  above  integral,  Eq.  (B3), 


becomes 


(B6) 


In  the  neighborhood  of  k*,  (3 (k)  ->  0(k*)  and  2u/h’(k)  assumes  an 
indeterminate  form.  Applying  L'HospitaPs  rule 

2u/h'(k)  2u7h”(k)  (B7) 

which  as  k  -4  k*  evaluates  to  (k  -  k*Vu  =  (2/h")1/2  and  this  result  in  turn 
implies  that  the  largest  acceptable  value  for  Ak  is  (2/h")1/2.  Thus,  the 
integral  in  Eq.  (B6)  reduces  to 


The  variable  of  integration  chosen  for  RAYMODE  is  obtained  by  setting 
y  -  (2/tc)1/2u,  so  that  Eq.  (B8)  becomes 


I  =  P(k*)/ti/a."(k*)  e~iIa<k  )+k  rJ{y*e’1y/2dy,  (B9) 

which,  of  course,  requires  that  a"(k)  does  not  become  vanishingly  small  as 

k  -4  k*. 

When  the  point  of  stationary  phase  falls  outside  the  segment  (ka,kb), 
the  phase  function  is  expanded  about  the  end-point  closest  to  k*.  In  this 
case,  the  Taylor  series  expansion  contains  a  first-derivative  term,  but 
through  a  suitable  transformation  the  integrand  is  cast  in  the  same  form 
as  the  one  in  Eq.  (B9).  Thus,  in  both  the  illuminated  case,  k*e  (ka,kb), 
and  in  the  shadow  case,  k*  e  (ka,kb),  the  integral  that  has  to  be 
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I 


I 


I 


I 


» 


» 


I 


I 


» 


evaluated  takes  the  form 


I 


-  fV< 


ig*/2 


dy. 


(BIO) 


This  integral  can  be  expressed  in  terms  of  the  Fresnel  integrals,  that  is, 

I  -  C(y2)  -  C(yi)  -  i[S(y2)  -  S(y1)]I  (B11) 


where 


and, 


v  2 

C(y)  =  Jocos|.x  dx  , 


(B 1 2) 


2 

S(y'<  -- j  sinS-x  dx. 


(B13) 


The  limits  of  integration  depend  upon  the  location  of  k*  relative  to  the 
interval  end-points  ka  and  kb.  Leibiger  (1971)  presents  a  table  of 
appropriate  values  of  yi  and  y2  under  various  conditions  on  ka  -  k*  and 
kb  -  k\ 


The  special  high  frequency  integration  (HFI)  method  is  essentially  an 
extension  of  the  original  RAYMODE  integration  (ORI)  method.  As  shown  in 
Sect.  45  of  Brekhovskikh  (1980),  when  a"(k*)  is  finite  the  integral,  Eq. 
(B1),  evaluates  to  the  "saddle  point  formula"  given  by 


I 


_  -fl°c0c*)  ±  ti/4] 

p(k*)  e 


(B14) 


where  the  sign  in  the  exponent  agrees  with  the  sign  of  a"(k*).  When 
|a"(k*)|  <  e,  for  some  small  e  >  0,  the  evaluation  point  (r,z),  say,  is  very 
close  to  a  caustic  -  an  indication  that  special  handling  is  required.  To 
cover  all  possibilities  the  phase  factor,  a(k)  +  kr,  is  expanded  to  an 
additional  term,  i.e., 


a(k)  +  kr  =  a(k*)  +  k*r  +  (a’  +r)(k  -  k*) 


+  a"(k  -  k*)2/2  +  a"’(k  -  k*)3/6,  (B15) 


) 


35 


where  the  derivatives  are  evaluated  at  k  -  k\  and  k*  is  a  point  of 
stationary  phase  that  lies  within  the  interval  of  integration  (ka,kb).  From 
the  stationary  phase  condition,  i.e.,  a'  +  r|k.  -  0,  the  a-derivatives  can  be 
expressed  in  terms  of  the  r-derivatives.  Let  rc  denote  the  value  of  range 
when  k  =  k*.  Then  at  k  «  k*.  a'  *  -rc,  a"  ■  -rc'  and  a"'  =  -rc".  Making  these 
substitutions  and  assuming  that  a"  =  0,  Eq.  (B15)  becomes 

a(k)  +  kr  s  a(k*)  +  k*r  +  (r  -  rc)(k  -  k*)  -  rc"(k  -  k*)3/6.  (B16) 

Under  these  conditions  the  integral,  Eq.  (B1),  can  be  put  into  the  form 


I  =  p(k*)  ©-ifct(k-)  ♦  k*r] 

x  J  expf  i[(r  -  rc)(k  -  k*)  +  rc"(k  -  k*)3/6]}.  (B17) 

Through  suitable  transformations  applied  to  the  variables  r  -  rc  and  k  - 
k*  (see,  for  example,  p.390  in  Brekhovskikh,  1980),  this  integral  can  be 
approximated  by  an  Airy  function. 

In  addition  to  the  extended  Taylor  series  expansion,  there  are  some 
subtle  differences  between  this  method  and  the  approach  taken  in  the  ORI 
method.  In  the  "original"  method,  the  derivative  r'  is  determined  using 
finite  differences,  whereas  in  the  "high  frequency"  method  it  is 
determined  analytically.  Finite  differencing  is  used  in  estimating  r" 
however. 
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Appendix  C 

4 

Breakdown  of  Phase  Factor 
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The  integration  procedure  identified  in  the  text  as  the  "Original 
RAYMODE  Integration"  method  is  applied  to  integrals  having  the  form  given 
in  Eq.  (A38),  i.e., 

kt>  -i[<x(k)  +  kr] 

I(kvkb)  =  j  (3(k)  e  dk.  (Cl) 


This  notation  conforms  with  that  used  by  Davis  and  Council  (1985b). 
Although  this  notation  is  convenient,  it  tends  to  disguise  a  feature  of 
some  importance  concerning  the  applicability  of  the  stationary  phase 
method.  Generally  speaking,  when  the  method  of  stationary  phase  is 
applied  to  integrals  of  the  form 

I  =fbA(x)eUf<X)dx.  (C2) 

0 


a  fundamental  consideration  pertains  to  the  amplitude  and  phase  factors 
A(x)  and  Af(x).  In  Eq.  (Cl)  the  amplitude  factor  is  a  rather  involved 
function  of  reflection  coefficient  amplitudes  and  WKB  amplitude 
functions.  Their  arguments  have  been  lumped  in  the  factor  a(k),  and  this 
factor  tends  to  "beat  against"  the  factor  kr.  Thus,  the  phase  factor  in  Eq. 
(Cl)  goes  through  positive  and  negative  swings  about  some  mean  value, 
whereas  the  amplitude  factor  -  even  though  fairly  complicated  --  varies 
smoothly  by  comparison.  Following  procedures  similar  to  those  that  lead 
to  Eq.  (B8),  the  integral  in  Eq.  (C2)  can  be  approximated  by 


A(x0)e 


1  Af(xa) 


c/*|f:i/2 


j 


e±1u2du 


-c/x|f"|/2 


(C3) 


for  any  small  e  >  0  (Bleistein  and  Handelsman,  1986).  This  approximation 
improves  as  e  0  and  eVA  — >  °°»  and  if  these  limits  are  approached 
Eq.  (C3)  can  be  replaced  by 


1(Af(x0)  ±  a/4] 

A(x0)  e 


(C4) 


► 
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where  x0  is  "the"  point  of  stationary  phase  and  the  assignment  of  +  or  - 
corresponds  to  the  sign  of  f"(x0).  Eq.  (C4)  is  referred  to  as  the  "stationary 
phase  formula."  The  analysis  that  leads  to  this  formula  is  referred  to  as 
the  "stationary  phase  method",  and  essentially  represents  a  particular 
approach  that  falls  under  the  more  general  heading  of  "saddle  point 
analysis"  (Bleistein  and  Handelsman,  1986).  A  tacit  assumption  made  ■> 
arriving  at  this  "formula"  is  that  there  is  only  one  point  of  stationary 
phase  falling  on  the  interval  of  integration. 

In  an  attempt  to  achieve  greater  clarity,  attention  is  focused  on  the 
a-component  of  the  phase  factor  in  Eq.  (Cl).  This  component  takes  the 
form 


where 


a(k)  -  y  +  2n  $(zu,Z|), 

V  -  4>(zuiz,)  ±  4>(zu,zs)  ±  <j>(zriz,). 


(C5) 

(C6) 


The  <j>  functions  are  of  primary  interest  here.  The  reader  is  reminded  that 
these  expressions  assume  the  relationship  zu  <  z9  £  zr  <  z,.  Reversing  the 
relationship  between  z8  and  zr  results  in  the  same  basic  expression  for  y 
but  with  the  roles  of  zs  and  zr  reversed.  The  <J>  functions  are  dependent  on 
mode  number.  Thus,  for  example,  the  full  expression  for  0(Zu,zs)  is 


£.5  —  - 

♦(zu,zp  m  =  fy/k  (z)  -  km  dz  . 

Zy 


(C7) 


If  each  wavenumber  is  multiplied  and  divided  by  some  reference  wave 
number,  say  k0,  this  integral  takes  the  form 


<l>m(Zu»Zs)  ■  k0C0  7m(ZuiZs)  ■  (0  Ym(zu>Zs). 


(C8) 


where 


1 

=  j 


C  J(z)  -  C"’  dz  . 


(C9) 


and  o)  -  27tf.  Note  that  Cm  is  the  "phase"  (or  ray  vertex)  velocity,  and 
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hence  using  Snell's  law  in  the  form  cos6  ■  C(z)/Cm.  this  integral  reduces 
to 


which  is  similar  to  what  Smith  (1974)  refers  to  as  a  "characteristic" 
time.  Indeed,  using  Smith's  terminology,  Ym(zu»zi)  's  one-half  Jm,  an 
"adiabatic  invariant"  for  mode  m.  Now  Jm  »  Tm  -  Rc(0mVCm>  where  Tm  is  the 
modal  cycle  time,  Cm  is  the  phase  (ray  vertex)  velocity,  and  co/Cm  -  km,  so 
that  the  exponential  (or  phase)  factor  can  be  expressed  in  the  form 

CI)[  (n  +  ±  Ym(zu,zs)  ±  Ym(zr>zl)  1 

+  kmr  -  (n  +  V2)  kmRc(0m).  (Cl  1 ) 

A  specific  assignment  of  the  ±  signs  can  be  associated  with  a  particular 
"ray"  path.  The  leading  term  in  brackets  [  ...  ]  does  not  oscillate  as  range 
increases.  The  remaining  terms,  however,  exhibit  significant  oscillations. 
These  terms  can  be  expressed  as 

km[  (k/km)r  -  (n  +  V2)  Rc(em)  ].  (Cl  2) 

Consider  the  first  cycle  (n  ■  0)  for  mode  m.  Then,  since  k/km  is  close 
to  one,  the  expression  in  brackets  (Eq.  (Cl  2))  is  approximately 

r-Rc(0m)/2,  (Cl  3) 

which  as  r  varies  from  zero  to  Rc(0m)  goes  from  a  value  that  is  negative  to 
one  that  is  positive.  Consider  the  simple  case  when  z8  -  zr.  The  mean  value 
of  phase  for  either  the  "up-up"  path  or  the  "down-down"  path  (i.e.,  Ym(zu*zs) 
and  Ym(zr«zi)  are  equal  and  opposite  in  sign  and  hence  cancel)  is  simply 
ojTm/2.  Thus,  the  phase  oscillates  about  coT m/2  with  swings  in  amplitude 
reaching  ±Rc(0m)/2.  Since  modal  cycle  distance  increases  with  increasing 


mode  number,  m,  these  oscillations  become  large  in  magnitude  for  high- 
order  modes  -  although  their  frequency  of  oscillation  becomes  less  rapid. 
Conversely,  osciMations  for  the  low-order  modes  are  not  as  extensive  in 
magnitude  but  they  are  fairly  rapid.  In  either  case  the  phase  variation  is 
rapid  relative  to  the  amplitude  variation,  and  hence  as  long  as  the 
reference  wave  number  k0  (  -  2itf/C0)  is  large,  i.e.,  the  frequency  is  "high", 
the  criteria  necessary  for  Eq.  (C4)  are  met.  The  validity  of  Eq.  (C3), 
however,  does  not  necessarily  require  that  k0  ->  <».  Even  for  relatively 
modest  values  of  k0,  as  long  as  the  interval  of  integration  is  properly 
confined,  this  integral  yields  an  acceptable  approximation.  Obviously  there 
is  some  frequency  below  which  this  method  fails.  Hopefully  the  transition 
to  the  RAYMODE  "normal  mode "  sum  occurs  before  that  frequency  is 
reached! 


Appendix  D 

Surface  Reflection  Loss  Submodel 
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The  surface  loss  submodel  used  in  RAYMODE  is  an  ad-hoc  adaptation  of 


results  derived  in  part  from  theoretical  considerations  and  in  part  from 
experimental  data  (Marsh  and  Schulkin,  1962).  That  is,  the  surface  loss, 
SL,  is  determined  from  an  expression  of  the  form 

SL-SL1+SL2,  (D1) 

where  SL,  is  given  by 

SL1  -  -20  log  (1  -  R)1'2,  (D2) 

and  where 

R  -  max  {  j  sine,  sine  -  (7tae2)'1/2exp(-ca02)sin8  }.  (D3) 

The  parameter  a  is  related  to  the  mean  square  slope  (jtan2p0)  through 

(2a)*1  -  jtan2p0  -  .003  +  .0026  Vs,  (D4) 

where  Vs  is  the  10-meter  wind  speed  in  kts  (Cox  and  Munk,1979).  The 
second  term  is  independent  of  angle  and  is  given  by 

Sl_2  -  -  10  log  {.3  +  .7/(1  +  .01(2x10-5  fVs2)2]}.  (D5) 


None  of  the  RAYMODE  documents  offers  a  derivation  of  either  of  these 
terms.  In  his  evaluation  of  RAYMODE  physics,  Deavenport  (1979,1982) 
notes  that  the  reflection  coefficient  can  be  determined  from 

R  -  1  -  \\  a  co$$r  d$r  d0r,  (D6) 

where  o  is  the  scattering  coefficient  evaluated  in  the  limit  of  large 
roughness.  This  expression  is  presumably  inspired  by  the  works  of 
Beckmann  and  Spizzichino  (1963)  and  Marsh  (1950).  Deavenport  then  notes 
that  a  is  proportional  to  the  mean-square  scattered  power  (as  given,  for 
example,  by  Eq.  (62),  p.89  in  Beckmann  and  Spizzichino.,  op.  cit.).  When 
the  indicated  integrations  are  performed,  the  final  form  of  SL!  is 
obtained. 


ft 


ft 
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Appendix  E 

Complex  Eigenvalues 


One©  the  appropriate  "reflection"  coefficients  have  been  determined, 
the  imaginary  component  of  the  eigenvalue  for  mode  m  can  be 
approximated  by  expanding  the  phase  integral,  evaluated  between  turning 
points,  about  the  real  part  of  km<  say  km",  i.e., 


dk 


m 


m 


where, 


Cy 

\J I<o/C(Z)l2  -  ki,  dZ  . 
Zu 


(El) 


(E2) 


Using  the  relationship  given  in  Eq.  (A34),  Eq.  (El)  can  be  expressed  as 

<j>(km)  =  0(km  )  -  V2  Rc(9m)  (^m  ’  km*)-  (E3) 

With  km  -  km*  +  iam,  then 

0(km)  2  '  iRc(®m)  am /2.  (E4) 

Therefore,  from  Eq.  (A32),  modal  attenuation  is  given  by,  approximately, 

am  2  ln|RuR,|/Rc(0m).  (E5) 

For  Z8  <  ZLand  Zr  <  ZL,  Ru  is  the  surface  reflection  coefficient  and  R,  is  the 
generalized  WKB  reflection  coefficient.  For  Z8  >  Zj_and  Zr  >  Z\_,  Ru  is  the 
generalized  WKB  reflection  coefficient  and  |R||  -  1. 


Appendix  F 

Departure  of  n2-linear  Profile  Model  from  Linearity 


Assuming  constant  salinity,  an  isothermal  BT  reading  is  usually 
interpreted  to  indicate  that  sound  velocity  varies  linearly  with  depth.  To 
see  that  the  n2-!inear  profile  model  is  generally  applicable,  consider  C(Z) 
for  Z  £  Zi_,  where  the  variation  in  depth  is  given  by, 

C(Z)  -  C0/(1  -  (FI) 

The  parameter  p!  can  be  determined  from  a  knowledge  of  the  surface 
sound  velocity,  C0,  and  the  profile  point  at  the  layer  depth,  (ZL,CL).  The 
velocity  gradient  dC/dZ,  from  Eq.  (FI),  is  given  by 


dC/dZ  -  i/2  p!  C3(Z)/C02, 

(F2) 

which  at  Z  -  ZL  is 

dC/dZ  -  i/2  Pi  Cl3/C02. 

(F3) 

Letting  AC  «  CL  -  C0, 

Pi  can  be  determined  from 

i-(c«/cL)2 

_  cl  -  cl  .  (2Cl+  AC)  AC 

(F4) 

1  zL 

Zi.C?  zLc  l 

and  hence 


dC  AC  (2  +  AC/CJ  9.  AC  ,  . 

dZ  2Zt  (1  -  AC/ZJ  Zt  ' 

which  is  just  the  definition  of  for  the  C-linear  (constant  gradient) 
profile  model.  Thus  at  Z  -  ZL  the  ratio,  say  p,  of  dC/dZ  for  n2-l inear 
variation  to  dC/dZ  for  C-linear  variation  can  be  expressed  as 

P  «  «VC.)  -  [£>££]’  *  C  -  (F6) 
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where  (3,  in  Eq.  (F3)  has  been  represented  as  2gi/C0.  For  C0  -  1500  m/s  and 
AC  «  1.8  m/s,  Eq.  (F6)  yields  p‘  -  1.0036,  which  suggests  that  the  n2-linear 
model  introduces  very  little  curvature. 

The  usual  demonstration  of  linearity  entails  expanding  the  radical  in 
the  denominator,  that  is, 

C(z)  =  C0[1  +  i/2  3iZ  +  3/a  pt2Z2  +  ...  ],  (F7) 

which,  since  «  Z,  reduces  to 

C(Z)  =  C0(1  +1/2p1Z)-C0  +  g1Z,  (F8) 

a  fairly  common  form  of  the  C-linear  profile  model. 


» 


> 
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Appendix  G 

Modal  Parameters  for  a  C3-linear  Duct 


I 


The  C3-linear  profile  model  has  a  depth  dependence  given  by 


C3 iZ)  -  C03(i  +  bZ),  Z  <  ZL 


(G1) 


where  b  -  3g0/C0  and  with  C0  and  g0  denoting  sound  velocity  and  it? 
gradient  at  Z  -  0.  The  WKB  phase  integral  for  a  surface  duct  is  given  by 
(see,  for  example,  Freehafer,  1951) 


Jc  ym  dz  -  (m  -  e)  2tc, 


(G2) 


where  ym  -  [k 2(z)  -  km2]1/2,  and  where  the  integration  is  carried  out  over  a 
full  cycle.  Applying  Snell’s  law,  this  integral  can  be  expressed  as 


fc  ym  dz  -  km  Jc  tane  dz. 


m  jc 


Since  cose  -  c/cm  -  (c0/cm)(1  +  bz)1/3,  then 


(G3) 


so  that 


dz  =  "  b^)  Sln0  COs20  d0' 

J  Vmdz  km  J  sin28  cos0d0 

c  c 


<G4) 


=  (4®)  stn30  ♦  const. 


(G5) 


The  integral  over  one  cycle  for  "mode"  m  is  equivalent  to  twice  the 
integral  from  0O  ("ray"  launch  angle)  to  zero  (where  "ray"  goes  horizontal). 
Assuming  that  the  highest-order  trapped  mode  has  turning  point  depth 
equal  to  the  layer  depth,  zL,  Eq.  (G2)  takes  the  form 


4nf  s*n0o  «  /M  i\„_ 
Jg-  — T~  =  ln 

rnc  ft. 


(G6) 


io  cos  0fl 
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1  L  L 

4  *  3  fl.  cos20. 


(G7) 


Since  dr/dz  -  cote  along  a  ray  trajectory,  then  the  cycle  distance  for 
"mode"  m,  say  Rc(0m),  can  be  determined  from 


v 

lRc(em)=  [  cote  dz 

2  V 


(G8) 


Substituting  Eq.  (G4)  yields 


jRc(em)  =  -|(^)  {  cos5ede. 


(G9) 


jV0™5  =  !(%■)*  s,n8« 0  *  isin20.) 

c?  Sine0  i  .2 

=  Ta*  ( 1  “  5Slne0). 


c*  9 

uo  ao 


(G 1 0) 


Let  rc  denote  cycle  distance  for  the  n2-linear  profile,  then  (e.g.,  see 
McGirr,  1983) 


jrc(0M>  =  T  sin0«cos0® 


(G 1  1 ) 


A  straight-forward  manipulation  casts  Eq.  (G10)  into  the  form 


(G12) 


For  ZL  -  50  m,  C0  -  1500  m/s  and  CL  -  1501  m/s,  the  quantity  in  square 
brackets  evaluates  to  1.00222  -  suggesting  that  the  two  profile  models 
produce  essentially  equivalent  "modal"  cycle  distances. 
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Appendix  H 

Sucker’s  Scattering  Integrals 


Bucker  (1980)  takes  surface  scattered  energy  into  account  by  adding 
randomly-phased  scattering  terms  to  the  coherent  mode  sum.  He 
essentially  uses  ray  theory  in  combination  with  mode  theory  to  couple 
surface-scattered  energy  into  the  total  sound  field  solution.  The 
additional  Mmode-to-ray"  and  "ray-to-mode"  interaction  terms  are 
incoherently  summed  to  get  the  total  intensity  due  to  the  scattered  field. 
The  component  of  intensity  due  to  scattering,  say  ls,  is  given  by 

U  =  51m  [Arin(Zs)Jm(Zr)  +  Am(Zr)Jm(Zs)]  6xp(-amr),  (HI) 
where  Am(Z)  -  |Urn(Z)/(km1/2Nm)|2,  am  -  2  lrn{km},  and 

Jm  *=  J[o(0m,0r)  exp(amp)/sin0r]  d0r,  (H2) 


and  where  p  «  r  -  rs.  The  reader  interested  in  specific  details  pertaining  to 
the  various  "modal"  quantities  appearing  in  these  expressions  should 
consult  the  original  article  by  Bucker  (ibid.,  1980). 


The  essential  step  toward  acquiring  an  appreciation  of  Bucker's 
approach  is  to  understand  the  processes  that  lead  to  "scattering"  integrals 
of  the  form  given  in  Eq.  (H2).  The  intensity  associated  with  the  mth  mode 
incident  on  the  surface  is  Im/cos0m.  Treating  the  elemental  scattering 
region  of  the  surface  as  a  point  source,  the  intensity  reduction  at  the 
receiver  is  given  by 


f-i_  cos0s/r 
r  |sin8dr/d0s| 


(H3) 


where  fr  is  what  Brekhovskikh  (1980)  refers  to  as  the  focusing  factor  for 
a  ray  launched  with  angle  0S  from  the  "source."  Let  o  be  the  scattering 
coefficient,  then  the  elemental  intensity  at  the  receiver  is 


d'sm  -o(lm/cos0m)  fr_1r8  dr8, 
or,  equating  0S  -  0m, 

cos8s/r 

r  |  sin0dr/d0s| 


(H4) 


(H5) 


Bucker  (1980)  expresses  lm  in  terms  of  modal  quantities,  i.e., 


Im  m  (2x/rs)  Am(20)  exp(-amrs),  (H6) 

and  then  integrates  Eq.  (H5)  to  get 

ism  “  (2rt/r)  Am(z0)Jm(z)exp(-amr),  (H7) 

where  Jm(z)  --  not  to  be  confused  with  a  Bessel  function  -  is  given  by 
Eq.  (H2).  The  limits  of  integration  are  defined  as  the  limiting  values  of  the 
set  of  rays  that  directly  connect  the  surface  and  either  the  source  or  the 
receiver.  That  is,  there  are  actually  two  interaction  (or  scattering) 
integrals.  Near  the  source,  rays  propagate  energy  to  the  rough  surface, 
whereupon  the  incoherent  component  of  the  surface  scattered  energy  gives 
rise  to  normal  modes  which  continue  the  propagation  of  waves  toward  the 
receiver.  This  interactive  effect  represents  a  ray-to-mode  exchange  of 
energy.  Near  the  receiver,  the  energy  associated  with  a  given  mode  is 
scattered  by  the  rough  surface,  whereupon  the  incoherent  component  of 
the  surface  scattered  energy  gives  rise  to  rays  that  propagate  energy  from 
the  surface  to  the  receiver.  This  interactive  effect  represents  a  mode-to- 
ray  exchange  of  energy.  The  details  of  this  procedure  are  omitted  here  but 
may  be  found  in  the  report  by  Gordon  and  Bucker  (1984). 

In  his  original  article,  Bucker  (1980)  assumes  that  the  specular 
surface  reflection  coefficient,  Sm,  is  given  by 

Sm  -  -exp(-gm/2).  (H8) 

The  parameter  gm,  referred  to  as  the  Rayleigh  roughness  parameter,  is 
defined  as  (e.g.,  see  p.  82  in  Beckmann  and  Spizzichino,  1963) 

9m  -  (2k«  h  sinesm)2,  (H9) 

where  ks  « ©/C8,  08m  is  the  surface  grazing  angle  of  incidence  for  the  mth 
mode,  and  h  is  the  rms  surface  roughness.  He  also  assumes  Lambert's  law 
(Houston,  1915)  which,  when  the  energy  lost  to  specular  reflection  is 
taken  into  account,  results  in  the  following  expression  for  o: 

a(0m.es)  -  V2  (1  -  ISml2)  sine8.  (H10) 
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